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Abstract 

The image reconstruction problem consists in finding an approximation of a 
function / starting from its Radon transform Rf. This problem arises in the 
ambit of medical imaging when one tries to reconstruct the internal structure 
of the body, starting from its X-ray tomography. The classical approach to 
this problem is based on the Back-Projection Formula. This formula gives 
an analytical inversion of the Radon transform, provided that all the values 
of Rf are known. In applications only a discrete set of values of Rf is given, 
thus, one can only obtain an approximation of /. Another class of methods, 
called ART, can be used to solve the reconstruction problem. Following the 
ideas contained in ART, we try to apply the Hermite-Birkhoff interpolation 
to the reconstruction problem. It turns out that, since the Radon transform 
of a kernel basis function can be infinity, a regularization technique is needed. 
The method we present here is then based on positive definite kernel functions 
and it is very flexible thanks to the possibility to choose different kernels and 
parameters. We study the behavior of the methods and compare them with 
classical algorithms. 



00 

in 



ii 



Contents 



Introduction and content v 

List of symbols vii 

1 Computed axial tomography 1 

1.1 X-rays 3 

2 Fourier based methods 7 

2.1 Characterization of lines in M 2 7 

2.2 The Radon transform 10 

2.2.1 Definition and basic properties 10 

2.2.2 Back projection 11 

2.3 The Filtered Back-Projection Formula 13 

2.3.1 The Central Slice Theorem 13 

2.3.2 The Filtered Back-Projection 15 

2.4 Filtering 16 

2.4.1 Filter resolution 17 

2.5 Discrete problem 19 

2.5.1 Phantoms 20 

2.5.2 Sampling 20 

2.5.3 Discrete filters 22 

2.5.4 Discrete functions 24 

2.6 Interpolation 28 

2.7 Discrete image reconstruction: Algorithms 29 

3 Algebraic Reconstruction Techniques 35 

3.1 Generation of the linear system 36 

3.2 Construction of A and p 38 

3.2.1 Algorithm 42 

3.3 Solving the system 43 

3.3.1 Kaczmarz's method 43 



iii 



iv CONTENTS 

4 Kernel based methods 47 

4.1 Hermite-Birkhoff interpolation 47 

4.1.1 Positive definite kernels 48 

4.2 Conditionally positive definite kernels 49 

4.2.1 Reconstruction by conditionally positive kernel functions 49 

4.3 Native function spaces 51 

4.3.1 Optimality of the interpolation method 53 

5 Kernel based image reconstruction 55 

5.1 Gaussian kernel reconstruction 56 

5.1.1 Regularization 58 

5.1.2 Regularization by truncation 60 

5.1.3 Regularization by Gaussian filtering 63 

5.2 Inverse multiquadrics reconstruction 65 

5.3 Multiquadrics reconstruction 67 

5.4 Compactly supported radial basis functions 69 

5.4.1 Compactly supported kernel reconstruction 71 

5.5 Scaled problem 75 

6 Numerical results 77 

6.1 Optimal parameters 77 

6.1.1 Window function parameters 78 

6.1.2 Kernel shape parameters 85 

6.1.3 Scale parameter 89 

6.2 Comparison of the methods 90 

6.3 Graphical user interface 92 

Conclusions 103 

Appendix A: Inverse multiquadrics kernel matrix 105 

Appendix B: Compactly supported kernel matrix 107 

Bibliography 111 



Introduction and content 



This thesis is the result of a three-months stage at the Univeritat Hamburg 
during which I studied the problem of clinical image reconstruction, i.e. the 
problem of obtaining the image of the internal structure of a sample start- 
ing from its X-ray tomography. From a mathematical point of view this 
correspond to find a function / knowing its Radon transform Rf. 

In the first Chapter the problem of image reconstruction is defined, we 
formalize the concept of Computed Axial Tomography and the history behind 
it. 

In Chapter 2 we discuss the mathematical aspect of the problem and its 
relation with the Radon transform. Then we follow the classical approach for 
solving the problem and deduce an inversion formula for the Radon trans- 
form: the Back- Projection Formula. Finally, we adapt the Back-Projection 
Formula to be used in real applications and thus we obtain the classical 
Fourier-based discrete image reconstruction algorithms. 

In Chapter 3 we introduce a different class of methods, called Algebraic 
Reconstruction Techniques (ART) and use them to solve our problem. 

Following the ART approach, in Chapter 4 we describe kernel based meth- 
ods and show how they can be used to solve the image reconstruction prob- 
lem. 

Chapters 5 and 6 are the original part of the work. In Chapter 5 we intro- 
duce a regularization technique that is necessary to implement kernel based 
image reconstruction and we realize such methods using specific positive def- 
inite kernel functions. In Chapter 6 we study from a numerical point of view 
the behavior of the methods in function of particular shape parameters and 
compare these methods with the Fourier based algorithms. 

In order to use the algorithms in a simple way, we also realized a graphical 
user interface that allows the user to test the algorithms on a set of predefined 
mathematical phantoms, with the possibility to choose options. 
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List of symbols 



l t> Q line in the plane characterized by values t and 6 

Rf(t, 6) Radon transform of the function / at a point (t, 6) 

S Schwartz space of rapidly decreasing functions 

Bh(x,y) back projection of the function h at point (x,y) 

F n f(u) n- dimensional Fourier transform of the function / at a point u 

F f(u) 1-dimensional Fourier transform of the function / at a point to 

R D f discrete Radon transform of the function / 

Boh discrete back projection of the function h 

F D f discrete Fourier transform of the function / 

FWHM(4>) full width half maximum of the function <f> 

\ y K(-,y) linear operator A applied to the function K with respect to the variable y 

R w f Radon transform of the function / multiplied by the window function w 

erf(x) error function evaluated at a point x 

Ff space of polynomial of degree lower or equal to k on R d 

k(A) 1-norm condition number of the matrix A 
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Chapter 1 

Computed axial tomography 



Computed axial tomography (CAT or CT) is a method that generates images 
of the interior of the body by digital computation applied to the measured 
transmission of X-rays tomography. In this process, an X-ray source and 
a set of aligned X-ray detectors are rotated around the patient (see Figure 
1.2(a)). The word tomography is derived from the Greek tomos (slice) and 
graphein (to write). 

The history of CT scan starts in Germany in 1895, when Wilhelm Conrad 
Rontgen (1859-1923; Figure 1.1(a)) discovered a new type of radiation, which 
he called X-rays [19]. This type of electromagnetic radiation, which has 
shorter wavelength then visible light and the ability to penetrate matter, 
was immediately used to image the interior of the human body. Figure 
1.1(b) shows one of the first X-ray images, this kind of images showed a two 
dimensional projection of the inner structures. In 1901 Rontgen received the 
first Nobel prize for physics. Basic to the CT technology are the theoretical 
principles of reconstruction of a three-dimensional object from multiple two- 
dimensional views relying on a mathematical model formulated by Johann 
Radon (1887-1956) in 1917 [17]. 

In 1979 the Nobel Prize for Medicine and Physiology was awarded jointly 
to Allan McLeod Cormack (1924-1998) and Godfrey Newbold Hounsfield 
(1919-2004), the two scientists primarily responsible for the development 
of computerized axial tomography in the 1960s and early 1970s. Cormack 
developed certain mathematical algorithms that could be used to create an 
image from X-ray data [2]. Working completely independently of Cormack 
and at about the same time, Hounsfield, a research scientist at EMI Central 
Research Laboratories in the United Kingdom, designed the first operational 
CT scanner, the first commercially available model and presented the first 
pictures of a patient's head [9]. Compared to a plan X-ray image, the CT 
image showed remarkable contrast between tissues with small differences in 
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(a) Picture of Wilhelm Rontgen (b) The first x-ray Frau Ront- 

gen's left hand 

Figure 1.1: The discovery of X-ray 



X-ray attenuation coefficient (Figure 1.2(b) shows the CT scan of a section of 
the brain). Since 1980, the number of CT scans performed every year in the 
United States has risen from about 3 million to over 67 million (for further 
details about X-ray history one should refer to [4] or [5]). 

The problem behind CT scans is essentially mathematical: if we know 
the values of the integral of two- or three- dimensional function along all pos- 
sible cross-sections, then how can we reconstruct the function itself? This 
is a particular case of what is called as an inverse problem and it was stud- 
ied by the Austrian mathematician Johann Radon in the early part of the 
twentieth century. Radon's work incorporated a sophisticated use of theory 
of transform and integral operators. 

The practical obstacles to implementing Radon's theories are several. 
First Radon's inversion methods assume knowledge of the behavior of the 
function along every cross-section, while in practice only a discrete set of 
cross-sections can be sampled. Thus it is possible to construct only an ap- 
proximation of the solution. Second, the computation power needed to pro- 
cess a multitude of discrete measurements and obtain from them a good 
approximation of the solution has been available for just a few decades. In 
order to overcome these obstacles theoretical approaches and approximation 
methods have been developed. 
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(a) A modern CT scanner (b) Brain CT scan 



Figure 1.2: Computed axial tomography today 

1.1 X-rays 

A CT scan is generated form a set of thousands of X-ray beams, consisting 
of 160 or more beams at each of 180 directions. When a single X-ray beam of 
known intensity passes through a medium, some of the energy present in the 
beam is absorbed by the medium and some passes through. The intensity of 
the beam as it emerges from the medium can be measured by a detector. The 
difference between the initial and final intensities tell us about the ability of 
the medium to absorb energy. 

The idea behind the CT scan is that, by measuring the changes in the 
intensity of X-ray beams passing through the medium in different directions 
and by comparing the measurements, we can determine which location within 
the sample are more or less absorbent than others. 

In our analysis of the X-rays behavior we will make some assumptions: 

• X-ray beam is monochromatic. That is each photon has the same en- 
ergy level E and the beam propagates at a constant frequency. If N{x) 
denotes the number of photons per second passing through a point x, 
then the intensity of the beam at the point x is 

I(x) = E ■ N(x); 

• X-ray beam has zero width; 

• X-ray beams are not subject to refraction or diffraction. 

Every substance has the property to absorbs a part of the photons that 
pass through it. To quantify this property we define the attenuation coeffi- 
cient of a material: 
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Definition 1. The attenuation coefficient of a substance is the fractional 
number of photons removed from a beam of radiation per unit thickness of 
material through which it is passing due to all absorption and scattering 
processes. 

In radiology a variant of the attenuation coefficient is used: the Hounsfield 
unit. Developed by Godfrey Hounfield, the Hounsfield unit represents a 
comparison of the attenuation coefficient of the medium with that of water. 
Specifically: 

Definition 2. The Hounsfield unit of a medium is 

jj _ ^medium ^water 

medium ~1 ; 

^iwater 

where A denotes the attenuation coefficient. 

Suppose now an X-ray beam passes through some medium located be- 
tween the position x and the position x + Ax. Suppose A(x) is the attenua- 
tion coefficient of the medium located there. Then the portion of all photons 
that will be absorbed in the interval [x,x + Ax] is p(x) = A(x)Ax. The 
number of photons absorbed per second by the medium is then p(x)N(x) = 
A(x)N(x)Ax. Multiplying both sides by the energy level E of each photon, 
we see that the loss of intensity of the X-ray over this interval is 

AI « -A{x)I(x)Ax. 

Let Ax — > to get the differential equation known as the Beer's law: 

£ = -A(x)I{x). (1.1) 

In other words: The rate of change or intensity per millimeter of a nonre- 
fractive, monochromatic, zero-width X-ray beam passing through a medium 
is jointly proportional to the intensity of the beam and to the attenuation 
coefficient of the medium. 

The differential equation (1.1) is separable. If the beam starts at the point 
x with initial intensity I = I(xo) and is detected, after passing through the 
medium, at the point X\ with final intensity Ii = I(xi), we get 



I ■ ill 

from which it follows that 



pi dl r xi 

/ — = - / A(x) dx, 



f Xl A(x)dx = In ) . (1.2) 
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Here we know the initial and final values of / and we want to determine the 
coefficient function A. Thus, form the measured intensity of the X-ray we 
are able to compute not the values of A itself, but the value of the integral 
of A along the line of the X-ray. 

From equation (1.2) it is easy to see that we can not discriminate two func- 
tions that have the same value of the integral along the X-ray path [x 0l x\]. 
The fundamental question of image reconstruction asks if it is possible to do 
that knowing the value of the integral of A along every line: 

The fundamental question of image reconstruction: Can we re- 
construct the function A(x,y,z) (within some finite region) if we know the 
average value of A along every line that passes through the region? (cfr. [4] 
pp. 7.) 

In our study of CT scans, we will consider a two dimensional slice of the 
sample, obtained as the intersection of the sample and some plane, which 
we will generally assume coincides with the xy-plane. In this context, we 
interpret the attenuation coefficient function as a function A(x, y) of two 
variables. 
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Chapter 2 

Fourier based methods 



In this chapter we study the methods that are used nowadays in the CT 
scanner. The mathematical foundation of these methods is based on the work 
of J. Radon on an integral transform, called in his honor Radon transform, 
and its inverse. Roughly speaking we can think that sending a set of X- 
ray beams through a sample and measuring the intensity of the beams after 
their passage through it, correspond to compute the Radon transform of the 
sample's attenuation coefficient. Thus applying an inversion formula of the 
Radon transform gives us the value of the attenuation coefficient within the 
sample. 

In theory this is possible if we know the value of the Radon transform 
in every point of the sample. In practice only a discrete set of values can 
be recorded by a X-ray machine, that's why we can only obtain an approx- 
imation of the original attenuation coefficient function and we will have to 
consider problems that arise working with discrete functions, such as sam- 
pling, filtering and interpolation. 

We start this chapter formalizing the concept of Radon transform. Since 
this operator involves the computation of the integral of a function along 
lines in the plane, we need first to define a suitable characterization of lines 
in R 2 . 

2.1 Characterization of lines in IR 2 

Consider again the equation (1.2): 

f 1 A{x)dx = In (^] . (2.1) 
Jxq \hj 

Suppose a sample of material occupies a finite region in space. At each point 
(x, y, z) within the sample, the material there has an attenuation coefficient 
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A(x,y,z). An X-ray beam passing through the sample follows a line I from 
an initial point P (assumed to be outside the region) to a final point Q 
(also assumed to be outside the region). The emission/detection machine 
measures the initial and final intensities of the beam at P and Q, from which 
the value ln(Jo//i) is calculated. According to (2.1) this is equal to the value 
of the integral Jpg A(x, y, z) ds, where ds represents arclength units along the 
segment PQ of the line /. Thus the measurement of each X-ray beam gives 
us information about the average value of A along the path of the beam and 
it is fundamental to find a useful representation of lines that can help us in 
solving the image reconstruction problem. 

For simplicity let assume that we are interested only in the cross-section 
of a sample that lies in the xy-plane. Each X-ray will follow a segment of a 
line in the plane and we look for a way of cataloging all such lines. 

The approach we adopt is characterizing every line in the plane by a 
point that the line passes through and a normal vector to the line. Then, let 
n be a vector that is normal to a given line /, then there exists some angle 9 
such that n is parallel to the line radiating out from the origin at an angle 9 
measured counterclockwise from the positive x-axis (Figure 2.1). This line is 
also perpendicular to / and thus intersects / at some point whose coordinates 
in the plane have the form (tcosO, tsinO) for some real number t. The line I 
is hence characterized by the values of t and 9 and so we denote / = l t fi. 




Figure 2.1: A line in the plane can be characterized by two real numbers i, 9 

Definition 3. For any real numbers t and 9, the line l t> g is the line pass- 
ing through the point (tcos^tsin^) and perpendicular to the vector n = 
(cos 6, sin 9). 

Because of the relationships lt,e+2-K = h,e and lt,e+-K = l-t,e for all t, 9, 
there is not a unique representation of the form l t $ for a line. For this reason 
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we will consider only the set of lines 

{kfi ■ t £ R, < 9 < tt}. 

If we consider the unit vector (— sin#, cos6>), perpendicular to n, every point 
on l t £ can be written as 

(t cos 9, t sin 9) + s(— sin 9, cos 9), 

for some number s £ K. So we can parametrize a line l t fi as (x(s),y(s)), 
where s£R and 

x(s) = t cos 9 — s sin 9 
y(s) = t sin 9 + s cos 9 

Note that for every point (x(s), y(s)) e Z t) g, we have x(s) 2 + y(s) 2 = t 2 + s 2 . 
With this parametrization the arclenght element along the line l t g is given 

by 

is) + (S) ds = ^(- sin9 ) 2+ ( cos9 y ds = ds 

Therefore for a given function A(x,y) dehned in the plane, we get 

/ A(x,y) =/ A(t cos 9- s sin 9,t sin 9 + s cos 9) ds. (2.2) 

The value of this integral is exactly what an X-ray emission/detection ma- 
chine measures when an X-ray is emitted along the line l t fi. 

Finally note that for an arbitrary point (x ,y ) in the plane and for a 
given value 9, there is a unique value of t such that (x ,y ) £ l tt g. The value 
of t is given by the solution of the system 

xq = t cos 9 — s sin 9 
y = t sin 9 + s cos 9, 

that is 

t = xq cos 9 + yo sin 9 
s = —Xq sin 9 + y cos 9. 

This formula will be used in the next sections to operate some change of 
variables that will be used in finding an inversion formula of the Radon 
transform. 
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2.2 The Radon transform 
2.2.1 Definition and basic properties 

The fundamental question of image reconstruction is: is it possible to recon- 
struct a function /, representing the attenuation coefficient of a cross section 
of a sample, starting from the value of the integral of / along every line l t fi 
in the plane? 

We will consider the integral of / for any values of t and 9, in other words, 
given a function / we associate to every point (t, 9) a number representing 
the value of the integral /j /. This leads us to the definition of the Radon 
transform: 

Definition 4. For a given function / : M 2 — > R, the Radon transform of / 
is defined by 

Rf(t,9) = / f ds — f(tcos0 — ssm$,tsm9 + scos9) ds, 
Ji t ,e 

Vt e R, 9 e [0,tt). 

So the Radon transform is an operator that, to a given function / of the 
Cartesian coordinates (x, y), associates a function Rf of the polar coordinates 
(t,9). 

Example 1. Consider a circle of radius r > and a function r r defined as 
follows: 

fr{ X ,y) = { 1 lf * 2 + y2 ^ 2 

( otherwise, 
since x 2 + y 2 = t 2 + s 2 , the Radon transform of / is 

Rf(t,9) — / f(t cos 9 — ssm9,tsm9 + scos^) ds = / Ids, 

JR J t 2 +S 2< r 2 

so we get 

vM-i'S^ .; f ; £ l "; 1 m 

[0 if t > |r|, 

Proposition 2.2.1. The Radon transform is a linear operator: for two func- 
tions f and g and constants a and f3, 

R(af + pb) = aRf + f3Rg. 

Proof. It follows from the linearity of the integral. □ 
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Example 2. Consider the function 



f{x,y) 



■ r 2 i 2^2 

if x + y < r 1 

•£ 2 - 2 i 2^2 

if r x < x + y < r 2 
otherwise, 



with < r\ < r 2 . We observe that / can be rewritten as / = f r2 — \ f ri , 
where f ri and f r2 are defined as in Example 1, by equation (2.3) and the 
linearity of the Radon transform, we get 



Rf(t,8) = Rf r2 (t,6)--Rf ri (t,9) 



' 2vV 2 2 - t 2 - sj- 

2^1 - t 2 




t 2 



if \t\ < n 

if ri < \t\ < r 2 
if \t\ > r 2 . 



Domain of Radon transform 

As we see from the definition, the Radon transform is a linear operator acting 
on functions and it involves improper integral on lines that can be infinity 
for some function. It is then natural to ask ourself for what kind of function 
is defined the Radon transform and in particular which is the domain of 
this operator, i.e. which space of functions is composed of all and only the 
functions that admit finite Radon transform. It can be proved (see [7]) that 
the space we are looking for is the Schwartz space 

S = if : sup | \x\ m P(d l7 d 2 )f(x)\ < oo, Vm e N, P polynomial 

of rapidly decreasing functions, but for the moment we can not consider this 
problem since the functions involved in medical imaging correspond to atten- 
uation coefficient of finite size samples and therefore are compact supported. 
In Chapter 4 we will face the problem of how to compute Radon transforms 
of functions that do not belong to S and we will find there some expedient 
to overcome this obstacle. 



2.2.2 Back projection 

Our aim is to recover a function /, representing the attenuation-coefficient 
of a sample, from the values of its Radon transform Rf. 

We start by considering a point (xo, yo) in the plane. For any values of 6 
there exists one and only one value of t such that the line l tt g passes through 
the point (x , yo)- I n particular, the value of t is t = x cos 9 + y sin 6. 
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In practice, any X-ray beam passing through a point (xo,yo) follows the 
line i( xo cos 0+j/o sin 0),0 f° r some angle 9. So the Radon transform Rf(xo cos 6 + 
yo sin^, 9) takes into account the value of the attenuation coefficient f(x ,y ). 

The first way one can try to recover f(xo,yo) is to compute the average 
of the Radon transform along all lines passing through (x ,yo)> that is 

1 f n 

— / Rf(x cos 9 + y sin 6, 9) d6 
7T Jo 

This leads us to the definition of the following transform, called back 
projection: 

Definition 5. Let h = h(t,9) a function in polar coordinates. The back 
projection of h at the point (x, y) is given by 

1 f" 

Bh(x,y) = — h(xcos9 + ysin9 1 9) d9 

7T JO 

Back projection is a linear transform: 

Proposition 2.2.2. The back projection is a linear transform, i.e. for all 
functions hi and h 2 and for all constants c\ and c-i, we have 

B(cihi + c 2 h 2 ) = CiBhi + c 2 Bh 2 

We observe that the back projection 

BRf(x, y) = - I' Rf(x cosO + y sin 9, 6) d,0 

7T Jo 

of the Radon transform, does not give us the value of f(x,y). Indeed, 
the value Rf(x cos 9 + ysin9,9) represents the total accumulation of the 
attenuation-coefficient / along a particular line. The integral BRf is com- 
puting the average values of those averages. Hence it gives us a smoothed 
version of /. 

Example 3. Consider fi a function corresponding to a disc of radius 1/2 
centered at the origin with constant density 1, that is 



1 if x 2 + y 2 < - 
otherwise. 



Then, for each line passing through the origin, we have Rf i(0, 9) = 1 and 
consequently BRfx(0,0) = 1. 
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Now suppose f 2 be defined by 




Again, for every line Zo,6» passing through the origin, we have Rf 2 (0,9) = 1 
and BRf 2 (0,0) = 1. 

Thus Bi2/i(0,0) = BRf 2 (0,0) = 1, but /i(0,0) = 1 and / 2 (0,0) = 0. 
This shows the fact that the back projection of the Radon transform does 
not necessarily reproduce the original function. 

2.3 The Filtered Back-Projection Formula 

In this section we will discuss the relationships between the Radon transform, 
the back projection and the Fourier transform. Thanks to these formulas we 
will obtain the inversion of the Radon transform. In other words we will 
be able to get the values of a function /, representing for example an X-ray 
attenuation coefficient, starting form the values of its Radon transform. 

In the next paragraphs we will consider successive transforms of a func- 
tion, for example in the central slice theorem in section 2.3.1 we will consider 
the Fourier transform of the Radon transform. In all these cases we will 
assume that all the transforms are well defined, i.e. we will assume that a 
function / belongs to the Schwartz space of rapidly decreasing functions S. 
For example one can think to / as a compact supported function. 



The interaction between the Radon transform and the Fourier transform is 
given by the Central Slice Theorem, also known as the Central Projection 
Theorem. 



We recall that the n-dimensional Fourier transform of a function / : 
R" R is defined as 



where i denotes the imaginary unit and x ■ uj the standard inner product in 
R™. For a function in polar coordinates f(t, 9), we consider the 1-dimensional 
Fourier transform F — F\, applied only to the variable t, i.e. 



2.3.1 The Central Slice Theorem 





We can now state the Central Slice Theorem in the case n = 2: 
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Theorem 2.3.1 (The Central Slice Theorem). For a function f defined in 
the plane and for all real numbers r, 9, 

F 2 f(r cos 9,r sin9) = F{Rf)(r,0). 

Proof. The definition of the Fourier transform gives 

/+oo r+oo 
/ f{x,y)e- ir( - xcose+ysind Uxdy (2.4) 
-oo J — oo 

Consider now the change of variables 

x — t cos — s sin 9 j t — x cos 9 + y sin 9 
y — t sin 9 -\- s cos 9 | s = — xs'm9 + ycos9. 

Note that the quantity t — x cos 9 + y sin is exactly the line / tj e . Moreover 
dxdy = dtds, indeed 

dx dx 

St ds 
oy oy 

dt 



cos 9 
sin 9 



— sin 9 
cos 9 



= cos 2 9 + sin 2 9 = 1. 



The integral in (2.4) becomes then 

/+oo r+oo 
/ f{t cos 9 - s sin 0, t sin 9 + s cos ^e^* ds dt = 
-oo J — oo 

/+oo / /-+oo \ 
I / f(t cos 9 — s sin 9, t sin 6* + s cos 9) ds ) e~ irt dt, 
-oo \J — oo J 

where we have factored out the inner integral since the term e~ irt does not 
depends on s. Now the inner integral in the last equation is exactly the 
definition of the Radon transform of the function / evaluated at point (£, 9). 
Thus the last integral equals 

/+oo 
Rf{t,9)e~ irt dt, 
-oo 

that is the definition of the 1-dimensional Fourier transform of Rf at the 
point (r, 9). 

In conclusion 

F 2 f{r cos 9,r sin9) = F(Rf)(r,9). 

□ 
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2.3.2 The Filtered Back-Projection 

Applying the back projection to the Radon transform gives a smoothed ver- 
sion of the original function. The following theorem, called Filtered Back- 
Projection Formula, shows how to correct the smoothing effect and recover 
the original function. 

Theorem 2.3.2 (The Filtered Back-Projection Formula). For all function 
f and for all real number x,y, 

f{x,y) = ^B{F~ 1 [\r\F(Rf(r, 9))]}(x, y). (2.5) 

Proof. By the Fourier inversion theorem, for any function / and any point 
in the plane (x, y), we have 

f(x,y) = F 2 - 1 F 2 f(x,y). 

Applying the definition we have 

1 /*+oo /*+oo 

f(x, y) = —J / F 2 f(X, Y)e^ x+ y^ dX dY. 

47T J — oo J — oo 

We pass now from Cartesian coordinates (X,Y) to polar coordinates (r, #), 
where X = r cos 9 and Y = r sin 6, with ret and 9 e [0, n] . Because of this 
change of coordinates in the integral, we have dXdY = \r\drd9 and so 

1 rTt z'+OO 

f( Xj y)= / F 2 /(rcos^,rsin^e ir(xcose+ysine) |r|^^- 

An 2 Jo J-oo 

Applying the central slice theorem to the factor F 2 f(r cos 9, r sin 9) = F(Rf(r, 9)), 
we get 

1 rir /*+oo 

f( X}V )= / F(Rf)(r,9)e ir{xcose+vs ' m(> ' ) \r\drd0. 

4tT 2 Jo J-oo 

In the last equation, the inner integral is by definition, 2tt times the in- 
verse Fourier transform of the function \r\F(Rf)(r, 9), evaluated at the point 
(xcos9 + ysin#, 9). So we can write 

f(x,y) = r F- 1 [\r\F(Rf)(r,9)](xcos9 + ysm9,9)d9, 
An Jo 

that is half of the back projection of the function F~ 1 [\r\F(Rf)(r, 9)]. Hence 
we finally obtain the desired formula 

f(x,y) = ^B{F- 1 [\r\F(Rf)(r,9)}}(x,y). 

□ 
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Observe that the factor \r\ in the formula (2.5) is fundamental. Indeed 
without this factor, the Fourier transform and its inverse, would cancel out 
and the result would be simply the back projection of the Radon transform 
of /, that as shown in example 3 does not lead to recover /. 

The Filtered Back-Projection formula is the basis for image reconstruc- 
tion. However it assumes that the values of Rf(t, 9) are known for all possible 
values (t, 9). In practice only a hnite number of X-ray samples are taken and 
we must approximate an image from the resulting data. 

2.4 Filtering 

Consider the Filtered Back-Projection formula in (2.5): 

f(x,y) = ^B{F- 1 [\r\F(Rf(r,e))]}(x,y) 

and suppose there exists a function (f>(t) such that F<fi(r) 
we could write 

\r\F(Rf)(r,9) = [F(j>.F(Rf)](r,9). 
By the properties of the Fourier transform we would have 

\r\F(Rf)(r,0) = F(<f>*Rf)(r,O), 

hence 

F-^rlFiRf^e)} = F- 1 [F( ( j ) *Rf)(r,8)} = 

= {<!>* Rf). 

We could then write equation (2.5) as 

f(x,y) = h( ( p*Rf)(x,y). (2.6) 

In this way the formula of the reconstruction of / would be simpler. The 
problem is that such a function <p does not exist. However, the previous 
discussion will be useful if we consider data Rf to be affected of noise, that 
is the case when we have to work with real data from the X-ray machine. 

Consider the function \r\F(Rf)(r, 9). The variable r represent a frequency 
that is present in a signal, so if the Radon transform has a component at 
high frequency, this component is magnified by the factor |r|. Since noise has 
high frequency, that means that the noise present in the image is amplified 
and this effect corrupt the reconstructed image. 



r . In this case 
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In order to obtain a formula less sensitive to noise, instead of \r\ we use a 
function, actually a low-pass filter, such that for r close to 0, it is near to the 
absolute-value function |rj, but vanishes if the value of \r\ is large. Moreover, 
in order to use the formula (2.6), in place of |rj we consider a function of the 
form A = F(f>, where A has compact support, or in other words, we consider 
4> band-limited function. 

In this way we obtain an approximation of /: 



Typically the function A is of the form A(co) = \uj\F(u)x[-l,l}{^>), for 
some L > 0, where %i represents the characteristic function of the set I. The 
function F is even and F(Q) = 1 in order to have an approximation of the 
function | • | near the origin and <p real valued. 

Typical low-pass filters used in medical imaging are: 

• The Ram- Lak filter: 



is simply the truncation of the absolute- value function to a finite inter- 
val. 

• The Shepp-Logan filter: 



f(x,y)^^B(F- 1 A*Rf)(x,y). 



A^uj) = \u\x[-l,l](w), 



A 3 (u) = \u\ y- 
( 2L 




sin(7rcc;/(2L))| if \uj\ < L 



= < 



7T 







otherwise. 



• The low-pass cosine filter: 



A 2 (uj) 



uj\ cos(ttuj/(2L))x[-l,l]- 



The plot of these filters in the case L = 10 is shown in Figure 2.2. 



2.4.1 Filter resolution 



Consider a function (j), suppose <fr > 0, with a single maximum value M in 
x = and increasing for x < 0, decreasing for x > (for example cj) can be a 
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Ram-Lak 

S hep p- Log an 




Figure 2.2: Main low pass filters used in medical imaging 

Gaussian). For another function /, the filtered version of / with (ft is given 
by / * (ft- 

Let now the numbers X\, x 2 be such that X\ < < x 2 and <j)(x\) = (ft(x 2 ) = 
M/2, half of the maximum value of (ft. The distance x 2 — X\ is called full width 
half maximum of the function (ft, in symbol FWHM ((ft). The resolution of 
the filter defined by the convolution with (ft is set to be equal to FWHM((ft). 




-10 -8 -6 4 -2 2 4 6 S 10 



Figure 2.3: Full width half maximum of a Gaussian 

To understand the reason of this definition, consider a function / that 
consists of two unit impulses separated by a distance of d. It is easy to 
show that if d > FWHM ((ft), then the graph of / * (j) has two peaks, but if 
d < FWHM ((f)), then the graph of f*(f> has only one peak and so we have lost 
of details. So we conclude that the smallest distance between two different 
features of / that can still be seen in the filtered signal / * <p is FWHM ((f). 
One should choose the filter function (f> in accordance with the resolution 
required. Intuitively we can think that a function (ft with a small FWHM is 
spiker than a function having large FWHM and has better resolution. The 
following examples can help to understand better 
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Example 4 (FWHM of a Gaussian). Let F(lo) = e~ B " , where wel and B 
is a positive constant. The maximum value of F is F(0) = 1. Half maximum 
is hence achieved for e~ Buj =1/2 i.e. u> = ±yJ\a{2)/B. therefore 



FWHM = 20n(2)/B 
Example 5 [FWHM of a the Lorentz signal) . The Lorentz signal is given by 

9[UJ) i + 47r 2 r 2 2 (cj- Wo ) 2 ' 

where u) 6 R and T 2 ,LOo are constants. The maximum of g is given by 
g(uio) = T 2 and g(uj) = T 2 /2 if and only if u — ±l/(27rT 2 ). Therefore 

FWHM = — 
vrT 2 




Figure 2.4: Gaussian filter (5 = 2, FWHM = 1.1774) and Lorenz signal 
(T 2 = 0.5, w = 0, FWHM = 0.6366). 



2.5 Discrete problem 

By Theorem 2.3.2 we know that if completed continuous data are available, 
then we can exactly reconstruct a function / starting from its Radon trans- 
form. In particular this is possible thanks to the back projection formula 
(2.5) 

f(x,y) = ^{F->|F(iZ/)(r,0)]}(x,y). 
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We have also seen, in Section 2.4, that in practice is convenient to replace 
the absolute- value function with a low-pass filter A. Thus, we may use the 
approximation 

f(x,y)*±B[F-\A*Rf)]{x,y). (2.7) 

In the practical implementation of this formula, we have to consider that 
only a finite number of values of Rf(r, 9) are measured by the X-ray machine. 
As a consequence of this fact we have to answer to some question about ac- 
curacy and computation. First of all we have to understand the sampling 
process, i.e. the process of computing only a discrete set of value of a contin- 
uous function; then we have to find the corresponding form of formula (2.7) 
for discrete functions; and finally we will use the process of interpolation to 
obtain value of the function we can not directly measure. 

2.5.1 Phantoms 

Different choices of filters, interpolation methods, and other parameters, will 
give us different reconstruction of the same image, thus we need a technique 
for testing the accuracy of one particular image reconstruction algorithm. 

In order to have a good accuracy test, we should know the original image 
we want to reconstruct. Moreover the method should be independent from 
the possible noise present in the data, but should depend only on the algo- 
rithm used in the reconstruction. To solve this problem, Shepp and Logan 
([21]) introduced the concept of mathematical phantom. 

A mathematical phantom (or simply a phantom) is a simulated object 
whose structure is defined by mathematical formulas. Thus no errors occur 
in collecting the data from the object and when an algorithm is applied to 
produce a reconstruction of the phantom, all inaccuracies are due to the 
algorithm. In this way we can compare different algorithms meaningfully. 

Figure 2.5 shows the well-known Shepp-Logan phantom. This phantom is 
widely used to test the quality of an image reconstruction algorithm since it 
is a good imitation of the human brain. 

2.5.2 Sampling 

Sampling is the process of computing the values of a function, or a signal 
defined in K, only on a discrete set of points {xk}kei- For example, points 
Xj~ can be taken with uniformly spacing, i.e. Xk = k ■ d for some positive 
number d, called the sampling spacing. The sampling spacing d determines 
the smallest detail of / that can be seen after sampling: if d is small we have a 
better resolution, while bigger values of d give us less resolution. On the other 
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Figure 2.5: The Shepp-Logan phantom 



hand, small values of d generate a bigger amount of data and make algorithm 
slower, so we want to find an optimal values of d which is a compromise for 
this trade-off. 

If we think to a signal as a sum of sinusoidal waves, the narrowest detail 
in the signal is given by the wave with the shortest wavelength (maximum 
frequency). If the signal is band limited, the Nyquist Theorem 2.5.1 below 
tells us that the signal can be completely recovered starting from its sampled 
version, provided that the sampling spacing is small enough. 

Suppose / band limited, i.e. its Fourier transform is zero outside a finite 
interval: F f(uj) = for \uj\ > L. If we extend Ff periodically out of [-L, L], 
its Fourier coefficients are given by 




thus 




J-L 



[ Ff(u)e lun i du = 2Lc_ 



Assuming Ff continuous we have 




Ff(uj) = c -» e 



—iumj^ 
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and so 



/(,)- *-■*■/<»)- £ /( ") ^-"*) 



that is / can be exactly reconstructed from the values f(nir/L), n 6 Z. 
The optimal sampling spacing is therefore d = 'j-, since L is the maximum 
value of \u\ in F/, the smallest wavelength is hence the optimal sampling 
distance is equal to half the size of the smallest detail present in the signal. 
This result is resumed in the following 

Theorem 2.5.1 (Nyquist Theorem). If f is a square integrable and band 
limited function, i.e. Ff(uj) = for all \oj\ > L, then for all x eK 

/(*)- t / (4) sil f (2.8) 



We observe that formula (2.8) involves an infinite series and that its 
general term sin (Lx — nir)/(Lx — nn) converges slowly. So we need a large 
number of samples for a good approximation. To address this, we can use a 
smaller sampling distance -|, R > L to gain a series with better convergence. 
This process is called over sampling. 



2.5.3 Discrete filters 

The image reconstruction formula (2.7) involves the inverse Fourier transform 
of a low pass filter. In practice also this function will be sampled like the 
Radon transform. Since the filters we consider are band limited, we use 
the Nyquist theorem 2.5.1 to know how many samples are needed to get an 
accurate representation of the filter. Here, we reconsider the filters introduced 
in section 2.4: 

• The Shepp-Logan filter is defined by 

, . . . /sin(7rw/(2L))\ , . 

' 2L 

— sin(7rcj/(2L)) if \uj\ < L 

< 7T 

otherwise. 

for some L > 0. The inverse Fourier transform of A s is a band limited 
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function and is given by 



1 f L 2L 

F~ l A i {x) = -l — sin (ttu)/(2L)) coscj duj = 

7T Jo IT 

L r/cos {Lx — it/2) cos (Lx + tt/2L)\ 
= tt 2 [\ x-tt/(2L) x + tt/(2L) ) 

1 1 
,x-?r/(2L) ~ x + tt/{2L), 

According to Nyquist theorem, F~ 1 A% can be reconstructed exactly 
from its values taken at distance ir/L. Setting x = mr/L, for n e Z, 
we get 

AT 2 

F " A ^ = 

The Ram-Lak filter is given by 

Proceeding as in the previous case we find that the inverse Fourier 
transform of the Ram-Lak filter satisfies 



1 



F~ l Ai{x) = 



7T 



Lxsm(Lx) 2sm 2 {Lx/2) 



x 2 



Setting again x = irn/L we obtain 

L 2 



F~ 1 Ai(irn/L) = 



2tt 



2sin(7rn) /sin(7rn/2)V 
irn \ nn/2 J 



Finally we consider the low-pass cosine filter: 

A 2 (uj) = \uj\cos(ttuj/(2L))x[-l,l}- 

The inverse Fourier transform of A 2 , evaluated at multiples of the 
Nyquist distance is 



F- l A 2 {nn/L) = 



2L 2 r7rcos(vrn) 2(1 + An 2 ) 
~^2~ [ I -An 2 ~ (1 - An 2 ) 2 



Figure 2.6 shows the sampled version of the inverse Fourier transform of these 
filters. 
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(a) Ram-Lak filter (b) Shepp-Logan filter (c) Cosine filter 



Figure 2.6: Sampled inverse Fourier transform 

2.5.4 Discrete functions 
Discrete convolution 

In order to implement formula (2.7) we have to decide what is convolution 
of discrete functions. 

A discrete function is a mapping from the integers into the set of real 
numbers. For a discrete function g, we write g n for g(n), for all n 6 Z. 

Definition 6. The discrete convolution of two discrete functions / and g is 
defined by 

+oo 

(/ *g)m= Yl fj9m-j Vm e Z. 
j=-oo 

The discrete convolution satisfies all principal properties of the standard 
convolution (e.g. commutativity and linearity). 

If only a finite set of values {/& = f(dk) : k — 0, . . . , TV — 1} is known, like 
in real applications, there exist two different ways to extend the sequence to 
all integers: 

1. Set f k = for all k t {0, . . . , N - 1}: 

2. Extend the sequence to be periodic with period N, f m = f m + n N, where 
for m 6 Z, n is the only integer such that m + nN e {0, . . . , N — 1}. 
We call such a function N -periodic discrete function. 

The convolution of two iV-periodic discrete functions is also a TV-periodic 
discrete function, defined by 

N-l 

{f *g)m= Yl fj9m-j Vm e Z. 
j=0 

Some problem can arise using discrete functions. For example if we are 
sampling a non periodic function, the periodic model is not the best to be 
used. But, even if the function is periodic, it may be not clear what the 
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appropriate period is and so we might sample the function on a set of values 
that do not correspond to one period, then extending data to form a discrete 
periodic function, we have the wrong one. The solution to these problems is 
a technique called zero padding. We take a finite set of values of a function g, 
then we pad the sequence with a lot of zeros and finally we form a periodic 
discrete function. 

The following theorem tells us that the convolution between a zero padded 
function and another discrete function gives the same result of " true" discrete 
convolution at least at the points where the values has been sampled. 

Theorem 2.5.2. Let f,g be discrete functions and suppose that 3K e N 
such that g k = for k < and k > K. Let M e Z, M > K - 1 and 
let f,g (2M + l)-periodic discrete functions defined by f m = f m , g m = g m 
for —M < m < M. Then for all m such that < m < K — 1 we have 
(f*g)m = (f*g) m - 

Remark 1. The proof of this theorem is just an application of the definition 
of convolution for discrete functions. For details we refer the reader to [4]. 

Discrete Radon transform 

In the context of a CT scan, the X-ray machine does not access the atten- 
uation coefficient along every line, but the Radon transform is sampled for 
finite number of angles 6 6 [0, n) and, for each angle, for a finite number of 
values of t. Values of 6 and t are equally spaced and we consider the parallel 
beam geometry: the X-ray machine rotates by a fixed angle and, at each 
angle, the beams form a set of parallel lines (Figure 2.7). 

If N is the number of angles at which the machine takes scans, then the 
values of 6 that occur are {k^, k — 1, . . . , N—l}. Assume that, at each angle, 
the set of parallel beams is composed of 2M + 1 equally spaced lines and let 
d be the distance between two lines, with the object to be scanned centered 
at the origin. Then the corresponding values of t are {jd : j = — M, . . . , M}. 

The continuous Radon transform Rf is then replaced by its discrete coun- 
terpart Rof, defined by 

Rohk = Rftid, kir/N) 

for j = -M, . . . , M and k = 0, . . . , N - 1. 

Theorem 2.5.2 above applies to the discrete convolution of the sampled 
band- limited function F~ 1 A and the sampled Radon transform Rof- Since 
the scanned object has finite size, we can set R D f(j, 9) = for \j\ sufficiently 
large. Thus with enough zero padding the discrete Radon transform can be 
extended to be periodic in the radial variable jd. 
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Figure 2.7: The parallel beam geometry 



For discrete function in polar coordinates, the discrete convolution is 
carried out in the radial variable only, so in the reconstruction formula (2.7) 
we have: 

N-l 

{F- l A*R D f) mfi = £ F- l A 3 R D f m _ ]fi . 

j=0 



Discrete Fourier transform 

Definition 7. The discrete Fourier transform of a TV-periodic discrete func- 
tion is another iV-periodic discrete function F D f defined by 

N-l 

(F D f)j = £ f k e- i2nki/N , for j = 0, . . . , N - 1 

and extended to be periodic for other values of j . 

The discrete inverse Fourier transform of / is given by 

(*S7)j = ^ E he i2nkj/N , for j = 0, . . . , N - 1 

fc=0 

and extended to be periodic for other values of j . 

The following theorems show that the properties of the Fourier transform 
are still valid for its discrete version. 



2.5. DISCRETE PROBLEM 



27 



Theorem 2.5.3. For a discrete function f with period N, 
FD 1 i F Df )n = fn, for all integers n. 

Theorem 2.5.4. For two N discrete functions f = {fk ■ < k < N — 1} 

and g = {gj. : < k < N — 1}, we have 

• F D {f*g) = {F D f){F D g); 

• F D (fg) = j f (F D f)*(F D g); 
. (F D f)j = (iW)-,; 

• Parceval equality: 

N-l i JV-1 

E \^\ 2 = ^ E \(Fnf)j\ 2 

3=0 iV j=0 

For a proof of these facts we suggest to see [4] . 
Discrete back projection 

In the continuous setting the back projection has been defined by 

1 

Bh(x,y) = — / h(xcos6 + y sin 6) d6. 

7T JO 

Now, in the discrete case, we replace the continuous variable 9 with angles 
kir/N for k = 0, . . . , TV — 1 and state the following 

Definition 8. The discrete back projection of a function h is defined by 

1 N ~ 1 7T 7T 

B D h(x,y) = — ^2 h(xcosk— + ysin/c— ,kir/N). 

In our case, B D has to be applied to h = (F^A) * (RdJ) and the re- 
construction grid within which the final image is to be presented is a rect- 
angular array of pixels located at (x m ,y n ), each of which is to be assigned 
a color or a gray-scale value. Hence B D needs the values of h at points 
(x m cos kn/N + y n smkir/N,k7r/N), while the Radon transform is sampled 
at points (jd, kit/N) arranged in a polar grid. The solution to this problem 
is interpolation. 
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2.6 Interpolation 

The process to obtain a function f(x), x 6 K starting form a discrete set 
of values f k = f(x k ), k — 1, . . . , iV + 1 is called interpolation. There exist 
several interpolation schemes. Here we give a short introduction to the most 
commonly used. 

• Nearest neighbor: f{x) = f k , where Xk is the closest point to x. This 
is the simplest method but generates a discontinuous function; 

• Linear: / is obtained connecting successive points (xk, fk), {x k +\, f k +\) 
with segment: 

f(x) = — —(x-x k ) + f k for x € [x k ,x k+1 ]; 
x k +i — x k 



Cubic polynomial spline: successive points (x k , f k ), (x k+i , f k+ i) are 
connected by apiece of a cubic polynomial. The pieces are joint together 
asking for C 2 continuity of the resulting curve. Also values of f'(x k ) 
are prescribed; 

Lagrange interpolation: / is given by a polynomial of degree N: 

j^l Uk^j (Xj - x k ) 



We notice that the nearest neighbor interpolation can be written as 
If(x) = £/"*[-!,!) (% ~m) , 

m u 

where \j denotes the characteristic function of a set J and f m is the value 
of the function / at the sample point md. Similarly the linear interpolation 
can be written If(x) = J2m / m A (| — m), with 



A(x) 



1 — |.x if \x\ < 1 
if \x\ > 1. 



Generalizing this approach we define, for a weighting function W satisfying 
certain conditions, the ^-interpolation Iw{f) of a discrete function / is 

M/) = E/^(i- m ) xeR - 
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We want I w f(kd) = f k . Then we choose W such that W(0) = 1 and 
W(m) = for all m e Z, m ^ 0. Moreover, if we want to preserve also the 
integral, we ask W to be such that 

/ Iw{f)(x)dx = d^fm- 

■JR m 

Then, W should satisfy 



f W(u) 



du = 1. 



Remark 2 (Interpolation and convolution) . Suppose that a discrete function 
g is given by the discrete convolution g — <f> * / and let W be a weighting 
function. Then the VF-interpolation 

W*) = m</> * /)(*) E E ^ - (^7^ - («* - *)) 

can be approximated as 

I w (<f> * f)(x) « Yl w (<f>)(x - kd)f(k), 
k 

that is, we can approximate the interpolation of (j) * / with a weighted sum 
of values f(k) and the interpolation Iw{4>) of </> at points x — kd (cfr. [4], 
pages 82-86). 



2.7 Discrete image reconstruction: Algorithms 

Having examined the discrete version of all elements in the formula (2.7), we 
have now all the necessary tools for approximating / starting from a discrete 
set of samples of its Radon transform. 

1. Image reconstruction algorithm I. Let I be the interpolation of 
(F^A) * (R D f), so that I(t, k-rr/N) is interpolated from the computed 
values (F^A) * (R D f)(jd, kir/N). Then for all points (x m ,y n ) in the 
grid, we approximate 

f(x m ,y n ) w ^B D I(x m ,y n ) = 

= —— / ( x m cos (k—') + y n sin (k—") , k—~) . 
2N ^ \ \ nJ y V Nr nJ 
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2. Image reconstruction algorithm II. Instead of interpolating the 
filtered Radon transform, we interpolate the filter and then, as shown 
in remark 2, we form a weighted sum of the sampled Radon transform: 

w ( k ) =Y, I f~ 1 a [x m cos (kjj) +y n sm (fc^J - jd,k^)R D f(jd,k^) 

j 

^ N-l 

f(x m ,y n ) -Wj^Yl w ( k )- 

We conclude this chapter applying the reconstruction formula in a par- 
ticular case. 



Example: crescent-shaped phantom 



We want to apply the reconstruction algorithm introduced in the previ- 
ous section to a particular phantom called crescent- shaped phantom (Figure 
2.8(a)) whose analytic expression is 



f(x,y) = < 



1 ifx 2 + y 2 <l A(x-i) 2 + y 2 >l 

2 lf 2 Y^ 
ifa: 2 + y 2 >-. 



In order to compute samples of the Radon transform, we calculate Rf ana- 
lytically. 

We observe that / can be written as a sum of two functions: f — fi — \f2, 
where fi and f 2 are given by 



fi(x,y) = 



1 if x 2 + y 2 < 
otherwise 



h{x,y) = 



1 if (x- -) 2 + y 2 < — 
v 8 ; J - 64 

otherwise 



for all (x, y) e R 2 . By the linearity of the Radon transform, we have 



Rf = Rh - - 2 Rh- 



(2.9) 



We know (see example 1) that for all fixed value r > 0, the Radon trans- 
form of the function 



fr(x,y) 



1 if x 2 + y 2 < r 2 
otherwise, 
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is given by 

and the function fi equals f ri for r\ = \. So we can use equation (2.10) to 
compute its Radon transform. 

Function f 2 is not of the form f r for some r, but can be obtained shifting 
such a function. More precisely f 2 (x,y) = f r2 (x — c,y), where r 2 — § and 
c = | . So we can use the shift property of the Radon transform for computing 

Rfx 

Theorem 2.7.1 (Shift property of the Radon transform). Let g : M 2 — > R a 

function and let G(t, 9) = Rg(t, 6) it's Radon transform. If 

h(x,y) = g(x - c x ,y - c y ), 

then the Radon transform H(t, 9) = Rh(t, 6) of h is given by 

H(t,9) = G(t - c x cos e-c y sin 9,9). (2.11) 

See [16] for more details about Radon transform shifting properties. 
Thus, from (2.10) and (2.11), we gain Rf 2 (t,9) = Rf r2 (t- c cos 9, 9), that 

is 

Rh{t e)= h\Jr 2 2 -{t-ccos9Y i{\t-ccos9\<r 2 
\ if \t - ccos^l > r 2 . 

By equation (2.9) we conclude that 



2yVf - t 2 if \t\ < ?~i A \t — ccos9\ > r 2 

Rf{t, 9) = <j 2 sfrl - t 2 - yfr\ -(t-c cos 9) 2 if \t - ccos9\ < r 2 
if \t\ > ri. 

Figure 2.8(b) shows the spectra of this function. 

Suppose now M = 20 and N = 18. We sample the domain [— 1, 1] x [0, 7r) 
with values tk = kd, k = —M, . . . , M and 9j = jj^, for j — 0, . . . , N — 1, 
where d = 0.05, obtaining the discrete Radon transform Rf\kd,jn/N). 

We consider Shepp-Logan filter as low pass filter and we consider ^ as 
sampling spacing. Indeed, when we used the continuous Fourier transform, 
we considered ^ as sampling spacing, in accordance with Nyquist theorem. 
Now, to compensate the additional factor 2n in the definition of the discrete 
inverse Fourier transform, we use ~ 2Z- To match this spacing with that 
of the Radon transform, we want ^ = d = 0.05 and so L = 10, then 



A(u) 



^|sin(0.057ra;)| |cj| < 10 
U > 10 
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(a) Phantom / (b) Radon transform Rf 

Figure 2.8: Crescent-shaped phantom and its Radon transform 



and 



(F -i A) ^ 



j (l-4n 2 )' 

Next we compute the discrete convolution 7 = F^ l A * Rf: for — 20 < 
rn < 20, < j < 17 

20 

-y(m,jir/N) = Yl (F^AU-kRf (0.05k, jn/N). 

k=-20 

Applying linear interpolation to the variable t of 7, we obtain 

20 

h(t,jir/N)= J2 l(m,jir/N)A(20t-m), t e [-1,1]. 

m=-20 

Finally we use the reconstruction algorithm I and we have the approxi- 
mation 



1 



17 



f(x,y) « — ^2 h(x cos (jir/N) +ysm(jir/N), jn/N). 



j=0 



Figure 2.9 shows the reconstructed function. 
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Figure 2.9: Reconstruction of / 



CHAPTER 2. FOURIER BASED METHODS 



Chapter 3 



Algebraic Reconstruction 
Techniques 

The Fourier based methods we have seen so far are the algorithms used in 
modern CT scan. Another approach to image reconstruction is based on 
linear algebra. Algorithms that use this approach are known as algebraic re- 
construction techniques, or ART. For example, the first CT scanner designed 
in the late 1960s by Godfrey Hounsfield used these methods. 

While the Fourier transform approach solves the continuous problem and 
then passes to the discrete one, ART considers the discrete problem from the 
beginning. 

Let us start reminding that an image is given by a grid of pixels (picture 
elements) and at each pixel is assigned a color (or a gray scale value) that 
represents the value of the attenuation coefficient in the region of the given 
pixel. 

Suppose that our image is formed by Kx K pixels, each of them represent- 
ing a small square in the plane. Define the pixel basis functions 61, 62, • • • , b K 2 

as 

f 1 if (x, y) lies in pixel number i 
bi(x, y) — { 

{ otherwise, 

and let Xi the color value of the i-th pixel. Then the resulting image can be 
written as 

K 2 

I(x,y) = ^bi(x : y)xi. 
i=i 

Applying the Radon transform to both sides, we get 

RI(t,9)=Y / Rb i {x,y)x l . 

t=i 



35 



36 CHAPTER 3. ALGEBRAIC RECONSTRUCTION TECHNIQUES 



The X-ray machine gives us the value of the attenuation coefficient func- 
tion / for some finite set of lines j = 1,...,J. Let us denote by 
Pj = Rf(tj,9j) these values. We want to approximate the attenuation co- 
efficient / with image /, so we set pj = RI(tj,9j) and rj.i = Rbi(tj,9j), for 
j = 1 , . . . , J and i = 1 , . . . , K 2 , and we ask that 

if 2 

Pj = Ys X i r i,i> j = l,...,J. (3.1) 

i=l 

Thus we obtain a system of J linear equations and K 2 unknowns. This 
system is very large but spare and typically overdetermined or underdeter- 
mined. We need then specific techniques for the solution of such a system. 
Before looking to these methods, let us see in detail how to generate the 
linear system (3.1). 



3.1 Generation of the linear system 

In this section we consider the problem of generate the linear system Ax = p, 
i.e. we want to compute A and p starting from the values of the Radon 
transform Rf of an attenuation coefficient function / obtained from a X-ray 
machine working with parallel beam geometry. 

We know the values Rf(tk, 0i) with tk = kd, k = —M, . . . , M and 
9i — Iff, j = 0, . . . , N. We want to compute A = (r jyi ) i = 1, . . . , K 2 j = 
1, . . . , J where K 2 is the dimension of the reconstructed gray-scale image 
I = {PXi} i=1 ^ K 2, with the components %i of the solution of the system 
representing the color of pixel PXi] J = (2M + 1)N is the number of sam- 
ples (tj,6j) on which Rf is measured and = Rbi(tj,9j) is the Radon 
transform of the i-th pixel-basis function 6j, computed at point (tj,6j), with 
bi defined by 

1 if {x,y)(EPX l 
if (x,y)^PX % - 
In order to solve this problem we assume that: 



h(x,y) 



1. The support of the function / and the samples points (tj, 9j) are con- 
tained in the unit square [—1,1] x [—1,1], this implies that d — ^; 

2. The reconstructed image / also lies in [—1,1] x [—1,1] and its center is 
at the origin (0, 0). If we consider J as a matrix I(r, s), r — 1, . . . , K, 
s = 1, . . . , K, whose components are the values Xi of pixels PX i: the 
center is the pixel of indexes r = |_^lpJ' s = Lh^J- We denote 
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3. The K 2 pixels in I are ordered as follows: 





4. Considering I as a function I : 



i.e. 



i(x,y) = Eik(x,y)xi, 



the Cartesian coordinates of pixels are PXi = [xi,Xi + \) x (y i+1 ,yi\. 
Thus, we are considering a top-down, left-right enumeration of vertexes, 
in accord with the matrix indexing. Note that PXi includes the top 
horizontal side and the left vertical side, but not the right and the 
bottom sides (see Figure 3.1), exception are the pixels in the last row 
and in the last column of I that include all sides. We identify a pixel 
with the coordinates of its top- left vertex (a^,^); 



(xi,yi 



(xi,yi-i 




xi+i,yO 



Xi+i,yi-i) 



Figure 3.1: Coordinates of a pixel 



5. X-ray beams has zero width. 

Using these assumptions, we find that pixel PXi, determinate by coordi- 
nates (Xi,yi), given by 



Xi i, 

c 



Vi = Hi. 

c 



In particular 

• PXi = [xi, x i+1 ) x (y i+1 ,yi\, x i+1 = x t + c _1 , y i+l = y { - c _1 ; 

• If i e {K, 2K, . . . , K 2 } => PX, = [ Xi , x i+l ] x j/J; 

. if i e {K(K- 1) + 1,...,K 2 } PX i = [x i ,x i + 1 )x[y i+1 ,y i ]. 
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3.2 Construction of A and p 

Assume that we know (tj,6j) and pixel coordinates PXi = x 
(yi + i,yi\. What we want to do now is to compute = Rbi(tj,6j). Let 

rn r 12 ■■■ r lK 2 
A-(r u )-\ 7 Ta - *f \-(A\A',...,A"). 
rji r. n ■■■ r JK 2 

and p = Rf(tj, OA. 

For a fixed i (E {1, . . . , K 2 } we define i l e M 2M+1 x the matrix such 
that ^(r, s) = Rbi(rd, s§), with r = —M, ... ,M, s = 0, . . . , N - 1, i.e. the 
columns of A 1 represent values of Rbi for a fixed value of 9: 

Rbi{-Md,0) ■■■ Rbi{-Md,(N - 

Rbi(Md,0) ■■■ Rk(Md, (N - l)f ) 
If R € M 2M+1 x R N is the matrix containing data Rf(t k , 9 t ) and if we set 

p = R(:) 

then we have that the i-th column of A is 

A 1 := A l (:), 

where the operator (:) indicates the analogous Matlab operator (see [8]). 

The problem can therefore be reduced to computation of the columns of 
A 1 for a fixed i. 

Let i € {1, ... , K 2 } and 9 e [0, ir) fixed. For t G R let r(t) = Rb^t, 9). 
We observe that since 6j = 1 inside pixel PXi and 6j = outside (and since 
we assume X-ray beams to have zero width), the value of r is the length of 
the intersection between line l tt g and PXi. Indeed: 

Rbi(t,9) — / bi(tcos9 — ssm9,tsin9 + scos9)ds = 
Jr 

= I ds = m(C), 

J{s€R:(t cos 6-s sin 8,t sin 6+s cos 6)ePXi} 

where m denotes the Lebesgue measure on M. and C = {s e 1 : tcos9 — 
ssin9 e [xi,x i+1 ), tsin9 + scos9 e (yi+i,j/i]} = / t ,e n PX^. 
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To determine for which values of t the line l t< e lies in PXi, we consider 
lines that pass through the vertexes of PJQ. Let 

P\ = (xi, Vi) Pi = {x u y i+1 ) P 3 = (x i+1 ,y i+ i) P 4 = {x i+1 ,yi) 

and let th be such that the line l th: g passes through point the Ph, h — 1, 2, 3, 4 
(Figure 3.2). Since (x ,y ) G l Xo c OS e+y l)S m0,e, we have 

ti = Xi cos 6 + yi sin 6 t 2 = cos 9 + sin 9 

ts = Xi+i cos 9 + yi + \ sin t 4 = cos 9 + yi sin 

Moreover, to determine the length of the intersection i tj g n PXj, we need 
to know the intersections between and the sides of PXi. Let 

-£'12 — h,e n -PiP P23 — ^,0 n P2-P3 P34 = ^,0 n -P3P4 -£-14 = ^,0 n P\Pa- 

Let us compute for example Ei 2 : 

k,e = {tcos9 — ssm9,tsin9 + scos9) = (x(s),y(s)) 
the line through P\P 2 is x = Xi, so we want x(s) = Xi, =>• 

t cos 9 — Xi t cos 9 — Xi t — Xi cos 9 

s = — — , y{s)=tsm9 + — — cos6> = — — . 

sm 9 sin 9 sin 9 

In a similar way we find E 23 , -E34, £41: 

/ t-XiCos^N (t-y i+ i sin 6> 

£12 = Xi, — a £23 = I 2 > 

V sin 6/ / \ cost/ 

_/ i-x i+ iCos6>\ (t-yism.0 

-£^34 — Xj_|_l, ; £/l 4 — 



sin 6* 7 V cos 9 



Of course this values are different in the limit case cos 9 = or sin0 = 0, 
i.e. for 6 = 0, tt/2. In these cases intersections between l tt g and PXi coincide 
with the vertexes Ph- 

Depending on the values of 6, the behavior of l tt g and l thJ g changes (see 
Figure 3.2). Therefore one should distinguish the cases 9 6 [0, |), 9 e 
[f , §), 6 G [§, §7r) and 6 G [|vr,7r). Consider for example 6 G [0, f ). In this 
case we have t 2 <ti < f 3 < t 4 , therefore 

if t < t 2 V t > t 4 l tfi n PXj = 
if t 2 < t < t 4 k,9 n PXi = AB 

where AB is given by 
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(a) 6 € [0, f ) (b) 6 6 [f , f ) 

Figure 3.2: Lines t fc for 9 <E [0, f ) 



1. AS = £ 23 £ 12 if t 2 < t < h; 

2. AB = E 23 E U ifh<t< t 3 ; 

3. AB = E U E 14 if t 3 < t < U; 

4. In the limit case i 6 {K, 2K, . . . , K 2 }, i.e. if we are considering a pixel 
in the last column of /, we have to account that side P3P4 6 PXi, so 

if Xi = x K 2 A 9 = A t = t A =$> r = AT3 = P 3 P A = c" 1 . 



We notice that in case 2. AB(t) = AB{t\) = AB(t 3 ) = const. Moreover for 
= 0, only cases 2. and 4. are possible, then we do not need to compute 
1/ sin 9. The determination of AB in the others 3 cases is similar. What 
remains to do now is to compute the length of AB. 

Let us start considering AB = E 12 E 23 (see Figure 3.3). The coordinates 
of the points are 

/ t-XiCos9\ ft-y i+1 sm9 \ 

En = Xi, — £, 23 = I - a , Vi+i 

V sm9 ) \ cos 9 ) 

We notice that the triangle E\ 2 P 2 E 23 is rectangle in P 2 and that, by defini- 
tion, 9 = P 2 Ei 2 E 23 . By Pythagoras theorem 



P2E23 
sin 9 
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(b) 0e[f,7r) 



Figure 3.3: Computation of AB 



thus 



P 2 E 23 =x Et . - xp 



t - y i+ i sin 9 



23 ~ r2 COS0 
t — y i+ i sin 9 — Xi sin 6> 



cos6> 



t-t 2 
cos 



We conclude that 



12-G>23 



* *' , for0e(O,^). 
sin cos 2 



We now consider = E 2S E U '- as stated before, for all t £ [£1,^3], 
G [0, f ), ^23^i 4 (t) = E^EuiU), hence 



EozE 



P1P2 



23- c '14 



COS COS 



We can reduce the number of cases if we order t^, ft = 1, 2, 3, 4 in increas- 
ing order: i min < t mm2 < t max2 < i max . Thus, we conclude that 



• for 9 ^ 0, 1 we have 



Rk(t,e) 



if ^min ^ t ^ ^min 2 
r 2 if t m in2 <t< t max 2 
. ^3 if ^max2 < * < t 

max? 



(3.2) 



where 
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• for 9 = 



c 



,-i 



if ^min 

< t < t, 

otherwise 



V (xi = x K 2 At — i max ) 



Rk{t,0) = 



max 



(3.3) 







• for 9 



2 



{ 



C 



-1 



if t min <t<t 
otherwise 



V (yi = Uk 2 A i = t min ) 



o 



(3.4) 



3.2.1 Algorithm 



We now summarize in a pseudo-algorithm the principal steps involved in the 
computation of A and p. Operations are indicated in Matlab language. 

1. Input: 

R: (2M+1)N matrix representing the Radon data 
K: dimension of the output image 

2. Initialization: 

A=zeros((2M+l)*N,K~{2}) ; 
c=floor((K+l)/2) ; 

3. Compute pixels coordinates: 

Y=repmat([0:K,K+l,l]); 
X=Y' ; 

x=X/c-l; y=-Y/c+l; 

4. Compute values th = th(xi, yi, 9), h = 1,2, 3, 4: 

T=zeros(K+l,K+l,N) ; 
for j=l:N, 

T(: , : , j)=x*cos(theta(j))+y*sin(theta(j)) ; 

end 

5. Compute A: 

• For all i=l:K*K extract sub-matrix Ti of T containing th values 
corresponding to pixel PXi] 
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• Sort Ti in increasing order; 

• for all j=l : N, for all k=l : 2M+1, compute A % (tk, Oj) using equations 
(3.2), (3.3), (3.4); 

• Fill i-th column of A: 
A(: ,i)=Atilde_i(:) ; 

6. Compute p: 
P=R(:); 

3.3 Solving the system 

In the previous section we saw how to compute matrix A and the r.h.s. p 
of the linear system generated from the algebraic approach to the image 
reconstruction problem. In this section we discuss other methods useful for 
the solution of this system. 

We start observing that the matrix A can be very large. In fact, every 
sampling of the Radon transform produces an equation, while at every pixel 
in the output image is associated an unknown. Moreover the system can be 
typically both underdetermined (more unknowns then equations) or overde- 
termined (more equations then unknowns). Another important property of 
the system Ax = p is that the matrix A is sparse. Indeed every particular 
line ltj,Bj passes through relatively few pixels in the grid. Thus most of the 
values Tjk are equal to zero. 

In order to solve the system we will use two different methods depending 
on whether the system is overdetermined or underdetermined. In the hrst 
case we will use the least square approximation, that means that the solution 
x will be given by x = argminj|^4y — b\\. In the second case we will use an 
iterative method called Kaczmarz 's method that will be discussed in the next 
paragraph. 

3.3.1 Kaczmarz's method 

The Kaczmarz's method [14] is an iterative procedure for approximating a 
solution of a linear system Ax = p, A 6 R mxn ) p g R m . If we denote the 
i-th row of A and Pi the i-th component of p, we can say that a vector x is 
solution of Ax = p if and only if 



r j • x = pi \/i — 1, . . . , m. 
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We also notice that the set = {x e R™ : r^ - x — Pi} is an affine subspace 
of R™. The idea of Kaczmarz's method is to project an initial approximated 
solution Xq on all these affine spaces, generating in this way a sequence of 
vectors, each of them satisfies one of the equations j-j • x — p^. 

Definition 9. Let L Pjr = {x e R n : r • x = p}, for r e R™ and p <E R, 
an affine space, let iieR". The affine projection of u on L p r is the vector 
u £ Lp^ r such that 

\\u — u\\2 — min — k|| 2 . 

Proposition 3.3.1. The affine projection u of a vector u in the affine space 
L Pt r is given by 

r ■ u — p 

u = u 7, — [775 — r. 

IMI2 

The Kaczmarz's method proceeds as following. From an initial guess it 
computes its affine projection on the first affine space. This projection is 
then projected on the next affine space in our list and so on until the last 
affine space. These operations consist of one iteration and the result of this 
iteration become the starting point of the next one. In detail the algorithm 
proceed as follow: 

f. Select £0; 

2. for k — 1, . . . , K max , x° k _ l = Xk-\ (where K maK is the maximum number 
of iteration allowed); 

3. for i — 1, . . . , m, 

4 

4. x k = x%_ 1 . 

The sequence xq, X\, x-i . . . generated by the method converges to a vector x 
that satisfies Ax = p (see Theorem 9.14 in [4] and references there). However 
the convergence can be slow and a lot of steps are needed to get a good 
approximation. Moreover if the system has no solution, like in many image 
reconstruction applications, then the behavior of the sequence is not clear 
and can be chaotic. 

In the field of medical imaging the size of the system can be a serious 
problem, but, as we know, the matrix A is also sparse. This means that 
when we compute x\ from x 1 ^ 1 , we only change the components of x]^ 1 that 
correspond to non zero entries of r^. So we can increase efficiency storing the 
location of these entries. 



xV 1 - 



Pi 



(3.5) 
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Another fact that is connected to the nature of the reconstruction problem 
is that adjacent X-ray beams transmitted along line lt t g, for similar values of 
t and 9, will intersect many of the same pixels, thus the corresponding affine 
spaces will be almost parallel. As a consequence, the convergence is slow and 
a lot of iteration is needed to reach a good approximate image. 

We conclude this section introducing a variation of the Kaczmarz's method 
that involves the introduction of a relaxation parameters in the formula (3.5). 
Let A^fc be such that < \^ < 2, then we replace formula (3.5) with 

i— 1 

T i - r *-l _ \. , ri ' X k ~P\ 
•^k ^k y H,k II ||2 ' i" 

I r i 1 1 

The parameter can accelerate the convergence of an indeterminate sys- 
tem. Note that if A^ = 2, then the vector x\ is just the reflection of x]^ 1 
across Li and there is no improvement in the proximity to a solution. That's 
why we consider < A^ < 2. 
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Chapter 4 

Kernel based methods 



In this chapter we present another approach for solving the image reconstruc- 
tion problem based on kernel functions. Reproducing kernels have already 
been used in image reconstruction [18]. Here we use a different approach. 

As usual our data are the discrete Radon transform of a function / : 
R 2 — > R {Rf(tj, 0j)}!? =1 , from which we want to find an approximation of 
the function /. 

The basic idea is to seek for the approximation s of / in a functions space 
S with finite dimension n, that is S — span{si, S2, ■ ■ ■ , s n }. Thus a function 
seS can be written as s = J2]=i c j s j f° r some Cj £ R. 

Then we ask that Rs coincides with Rf on the points {tj,9j) for all 
j = l,...,n, i.e. 

{Rs)(t j ,0 j ) = {Rf)(t j ,e j ) j = l,...,n. (4.1) 

By linearity of R the coefficients Cj, j — 1, . . . , n are given by the solution of 
the linear system Ac = b, c 6 R™, where 

A = {Rs k {tj,9j))j,k=l,...,n b = {Rf(tj, 9j))j=l,..., n - 

4.1 Hermit e-Birkhoff interpolation 

We generalize the image reconstruction problem (4.1) considering the prob- 
lem of finding a function s 6 S such that /|a = s|a for some function /, 
where A = {Ai, . . . , X n } is a set of linearly independent linear functionals 
(see [12, 13]). In our specific case we will consider Xjf = Rf(tj,9j). We also 
assume S = span{si,S2, ■ ■ ■ ,s n } with |A| = n. By linearity, the problem is 
equivalent to the linear system 

Ac = / A , 
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where 

n 

A — (XjSk)j,k=l,...,n A = (Aj/)j=l,..., n S = ^ CfcSfc. 

j=l 

Theorem 4.1.1 (Mairhuber-Curtis [15], [3]). Let Sia d ,d>2, suppose 
$1 contains a interior point, then there is no Chebjichev system S\, . . . , s n , 
n > 2 on Q, i.e. for all s\,...,s n real valued functions on Q,, exists A = 
{Ai, . . . , A„} ; |A| = n such that the matrix (\k{sj))j,k=i,...,n * s singular, where 
Afe = for pairwise distinct points £ Q are the functionals "evaluation 
attk". 

This theorem tells us that if we want to find a basis si,...,s n of S such 
that the system (4.1) has a unique solution for all data A, the basis should 
depends on the location of the data, i.e. on A itself. For this reason we will 
choose 

Sj = \MK(;y), 

where K : R d x R d — > K and AJ indicates that operator Xj is applied to 
variable y. 



4.1.1 Positive definite kernels 

The problem is then solving /|a = s|a with s = J2j c jX V jK(-,y). It can 
be also written as a linear system A K ^c = /a, where c e R™ and A K a = 
(XtX]K(x,y))l k=1 . 

In order to have a unique solution for every choice of A, Ak,a must be 
non-singular. This is certain true if we assume K symmetric, i.e. K{x, y) = 
K(y,x) for all x,y 6 M d , and positive definite, that means that Ak,a is 
positive definite for all A. The following theorem gives us a characterization 
of positive definite functions. 

Theorem 4.1.2 (Bochner). Assume : R d — > M even and continuous and 
that its Fourier transform <I> is such that the Fourier inversion theorem holds: 



( Z7T ) JR d 



(2tt)° 

then, if <&(uj) > 0, Vcj e R d , K(x,y) = <J>(||x — y\\) is positive definite. 

Example 6 (Gaussian). <E>(x) = e - "^' 2 is positive definite since = 

_imi 2 _ 
e ~r- > 0; 

Example 7 (Inverse multiquadric) . $(x) = ^/=jj p i s positive definite since 

$(a>) = ifd-i (||a;||)||a;||~^ > 0, where if„ denotes the Bessel function of 
2nd kind of order v. 
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4.2 Conditionally positive definite kernels 

Definition 10. A set of functional A is said to be /c-unisolvent (or unisolvent 
w.r.t P^) if for p e Ff we have 

p\ A = p = 0, 

where Ff denotes the set of all polynomial of degreed less or equal of k on 
R d . 

Definition 11. The radial kernel $ is conditionally positive definite of order 
k and we write $ e cdp(k), if for K(x, y) = — y\\) the quadratic form 

n 

c t a ka c= x: Cic^mix-vw) 

is positive for all possible A, |A| = n and vectors cel" \ {0} satisfying 

n 

4.2.1 Reconstruction by conditionally positive kernel 
functions 

We introduce a polynomial part in the interpolant s, so what we have to do 
now is to solve /[a = s|a with s of the form 

n 

- = E^(||--y||)+p, 

where p e P^-i (A; = fc($)) and vector c £ M™ satisfying the vanishing 
moment condition 

n 

E^A J (p) = VpePjU 
J'=l 

Theorem 4.2.1 (Michelli,Wu). TTie reconstruction problem /|a = s|a 
under vanishing moment condition a unique solution s, provided that $ e 
cdp(k) and the Junctionals A are (fc — 1) -unisolvent. 

Proof. Let pi,...,p m £ Pfc_i a basis of Pfc_ l5 where m = dim(P^_ 1 ) = 

n m 

a = E^Ay*(ii--y||) + x;*p». 

i=i i=i 
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for some c £ 1" d 6 R m . Condition /a = s|a an d the vanishing moment 
condition, are equivalent to the linear system 

A?s = A?/ Vi = l,...,n 

m 

£ 9 A*(p z ) = W = l,...,m 
I j=i 



that in matricial form is 



* P \( c \ _ ( /Ia 

p t o U r o 



(4.2) 



where = AfAj$(||x - y||), i,j = l,...,n, P^, = A^p,), j = l,...,n, 
I = 1, . . . ,m and O e M m x M m and e M m are a matrix and a vector with 
all components equal to zero. 

We consider the homogeneous system 

$c + Pd = J c r $c + c r Pd = 
P T c = ^ | c T P = 

substituting the second equation in the first one, we have c T §c = 0. Since 
$ 6 cdp(k), c = 0. The hrst equation becomes then Pd = 0, but A is (A; — 1)- 
unisolvent, that implies d — 0. So the unique solution to the homogeneous 
system is the null solution. □ 

Example 8. Conditionally positive functions: 
f . Polyharmonic splines: 



<p(r) 



r 2k d log r if d is even 
r 2k ~ d if d is odd 



2k > d; 

2. Gaussian: tp(r) = e~ r2 , k = 0; 

3. Multiquadrics: p(r) = (1 + r 2 f , // > 0, i/ ^ N, k = [z/|; 

4. Inverse multiquadrics: </?(r) = (I + r 2 )", v < 0, A; = 0; 

5. Power function: <p(r) = r?, < (3 $ 2N, fc = [§]. 
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4.3 Native function spaces 

In this section we will show that the solution of the Hermite-Birkhoff inter- 
polation is optimal in the sense that it is the function of minimum norm 
among all functions that interpolate data /|a, where the norm is taken in a 
suitable Hilbert space. 

For a fixed positive definite function K : M. d x M. d — > M., we define the 
function spaces 

S A = span{\yK(., y) : A e A} S = {s € S A : | A| < oo} 

and the dual space 

n 

L = {\ = A CjA = cjXj : ceK", | A| < oo}. 

We observe that, for all s £ S there exists A e L such that s = s\ = \ y K(-, y). 
Indeed 

n n / n \ ^ 

s = £ CjX'jKi; y) = £ CjX'Ki; y)=[ y £ *A y) = A"if (-, y), 

3=1 3=1 \j=l J 

with A = J2 c j^j £ L. 

We define an inner product on L: 

n n 

(A, n) K = X x ^K(x, y) = YH Cjd k \^lK(x, y), 

j=l k=l 

where A = J2 c j^j an d /-t — £ ^fc/Wfe, and the norm 

||A|U = (A,A)f . 

Thanks to the duality relation between L and S 1 , we introduce a topology 
also on S so that (sa,s m )k = (A,^)at, ||-|U = (v)* 
Remark 3. L = S: L and 5 are isometric with respect to the norm 
Moreover, for all \i 6 L, /x is continuous on 5. In fact 

|/i(s A )| = (//^AT^y) = |(//,A)^| < ||/j|U||A|U = ||//|UIMU- 

We now set D = L and F = S the topological closures with respect to 
IHIif. The following theorem holds: 

Theorem 4.3.1 (Madych-Nelson, 1983). For all s M € S, for all A e D we 

have 

{\*K{.,y), Sll ) K = (s x ,s^) K = (\,fj) K = \ x » y K(x,y) = A(s„). (4.3) 
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Proof. The statements holds for A, [i e L and the representation (4.3) follows 
by continuity. □ 

Definition 12. Let H = {/ : S] C l" 1 -> R} a Hilbert space, K : tl x ft R 
is a reproducing kernel for i/ if: 

1. e # for all a; € ft; 

2. /(x) = (/, K{; x)) H , for all f £ H, x € ft. 

Corollary 4.3.2. K : R d x R d — > K is t/ie reproducing kernel of the Hilbert 
space F. 

Proof. We prove properties 1. and 2. of Definition 12 

1. For 5 Z e L, z£ R d , 5 y z K(-, y) = K(-, z) e F, for all z e R d ; 

2. For A = 5 Z e L, by theorem 4.3.1 (#(-,*),/)* = /(z), for all / e 

F, 2 e R d 

□ 

Corollary 4.3.3. The point evaluation 5 Z : F — > R are continuous on F. 
Proof. \S z (f)\ = < \\5 z \\ K \\f\\ K for all zj. □ 

Corollary 4.3.4. Let s = Sf t \ £ S denote the unique interpolation to f £ F 
on A: s\\ = f\\, then the Pythagoras theorem 

\\f\\ 2 K = \\s\\K+\\f-s\\ 2 K 

holds. 

Proof. For s = \ y K(-, y) £ S, A £ L, we find (s, g)x = for all g such that 
X(g) = 0, i.e. s = S/ ; a is orthogonal to the kernel of the functional A £ A. 
Hence 

(■Sf,A,f ~ S fA ) K = 0, 

i.e. the interpolant Sf,\ £ S is the orthogonal projection of / £ F onto S. □ 
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4.3.1 Optimality of the interpolation method 

The following results are consequences of Theorem 4.3.1 and its corollaries. 

Theorem 4.3.5. The interpolant s = Sf\ 6 S is the unique minimizer of 
the energy functional \\-\\k among all interpolants to data f\\, i.e. 

\\s\\k < \\g\\K Vg e S s.t. g\ A = f\ A . 

In this sense the interpolation scheme is optimal. 

Corollary 4.3.6. The interpolant s = s/a 6 S is the unique best approxi- 
mation to f £ F from S with respect to \\-\\k- 
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Chapter 5 



Kernel based image 
reconstruction 

In this chapter we apply the kernel based methods saw in Chapter 4 to the 
problem of image reconstruction. We will see that the Hermite-Birkhoff in- 
terpolation can not be applied to the original reconstruction problem because 
the Radon transform of a kernel basis function can be infinity. We will then 
overcome this obstacle introducing a regularization of the integrals involved 
in the computation of the Radon transform. Thanks to this regularization it 
is possible to generate a liner system, solving the linear system one can find 
an approximation of the image to reconstruct. 

This technique can be used with both parallel beam geometry and scat- 
tered data. This second case is useful when one wants to reduce the dosage of 
X-rays passing through the sample. Scattered data can then be interpolated 
using suitable methods. For example a radial functions method was recently 
introduced by Beatson and zu Castell [1] to obtain new values of the Radon 
transform. This technique can be combined with the methods introduced 
below to obtain a reconstruction using less initial data. 

Let / : R 2 — > R be a function. Consider again the problem s\\ = /|a, 
where 

f\ A = {Xjf}] =1 , Xjf = R[f(x 1 ,x 2 )](t j ,6 j ), j = 1, . . . ,n, 

and s is an approximation of / belonging to the space 

S = ^2cjX^K(-,y) : c 3 E R, K(x,y) = 3>(||x - y\\) positive definite J . 

Let us denote b 3 {x) = X y j (K(x,y)) = R[K(x,y)](t j ,e j ), x e R 2 , 1 < j < n, 
the basis of S, so that s(x) = YIj=i c jbj{%) for some c 6 R™. The interpolation 
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conditions s\a = /|a a re equivalent to X k s = X k f Vfc = 1, . . . , n. By linearity 
of the Radon transform we obtain 

n 

Yl c 3 x m K ( x ^ y) = x kf, k = i,...,n, 

or in matrix form 

Ac = f A (5.1) 

with f A = (Ai/, . . . , X n fY and A = (ajy)i<fcj<„ given by 

a kJ = X%X v jK(x,y) = R x {R y [K(x,y)](tj,9j)} (t k ,9 k ). 

Thus, to determine c and then solution s, we have to solve the linear system 
(5.1). 

The first step in solving system (5.1) is of course computing the matrix A. 
We start by considering a generic basis function bj(x) = R y [K(x,y)](tj,9j) 
(for simplicity of notation we omit index j and so we denote (tj, 9j) = (t, 6)). 

We notice that since the kernel K is of the form K(x, y) = Q(\\x — y||), we 
can use the shift property of the Radon transform to simplify the computation 
of bj. Indeed, if we set k(y) = K(0,y) = $(\\y\\), then K{x,y) = <E>(||x - 
y||) = k(x — y) = k(y — x). Hence, by theorem 2.7.1 (shift property), if 
g{t,e) = R[k(y)}(y,9), we have 

R}>[K(x,y)}(t,6) = R»[k(y-x)](t,0) =g(t-x-v,0), 

where v = (cos 6, sin 9). So, in order to obtain bj we have only to compute 
R[k(y)](t,9) = R[K(0, y)](t, 9). Notice that this property is independent of 
the particular kind of kernel (Gaussian, multiquadrics, etc.) used and so is 
applicable with any kernel function of the form K(x, y) = $>(\\x — 2/11)- 

5.1 Gaussian kernel reconstruction 

We start considering the Gaussian kernel 

K(x,y) = e -M 2 . 

R[k(y)](t, 9) — j k(t cos 9 — s sin#, t sin 9 + scos9) ds = 




exp (—(9 — s sin(2) 2 — (t sin 6* + s cos9) 2 ) ds = 



= [ e -(' 2 +* 2 ) ds = e"* 2 / e" s2 ds = 
= v^e"' 2 =g(t,9). 
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Thus 



Ry[K(x, y)](t, 9) = g(t-x- v, 9) = v^e"^ 2 
and we conclude that 



bj(x) = \fv:e 



-(tj—X-Vj) 2 



where Vj = (cos 9j , sin 6 3 , 



We now want to compute a kj = R[bj](t k ,9 k ). Again for simplicity of 
notation, we write (tj, 9j) = (t, 9) and (t k , 9 k ) = (r, ip), then 

R[s/Tre~ < ' t ~ x ""^ 2 ] (r, (p) = s/n / exp (— [t — (r cos ip — s sin (p) cos 9 

Jr 

— (r sin p + s cos ip) sin 9} 2 ) ds = 

= -\/tt / exp (— [t — r(cos <p cos 9 + sin ip sin 9)+ 

Jr 

+s(sin ip cos 9 — cos (/? sin 6 1 )] 2 ) ds = 

= / exp (—[t — r cos (9? — 9) + s sin (<p — 9)] 2 ) ds. 
Jr 

If we set a = sin (ip — 9) and b = t — r cos ((p — 9), we can write 

a fc , = J R[ v ^e-( t - a; - ?, ) 2 ](r,^) = ^ [ e^ as+b ^ ds. 

Jr 

Hence, if a 7^ we have 



a kj = ^ [ e -(-+^) 2 d(as + b) = ^ / e"" 2 du = -, 
a </r a a 

while, in the case in which a = 0, 

•/1 

since both ip and are in [0, it), a = if and only if 9? — 9, so we conclude 



a k,j — \ 



7T 



+ 00 if ^fc = % • 
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5.1.1 Regularization 

We saw that matrix A can have infinity entries. More precisely R[bj] (t k , 6k) = 
+oo for some values of j and k (those values s.t. 6k = 6j), that means that 
for these values bj(x) is not integrable on line k k ,e k - 

To overcome this obstacle we must find some regularization technique so 
that the value of the Radon transform of the basis elements bj is finite for 
all k,j. The simplest choice is to consider a truncation of the integral, i.e. 
computing 



instead of the integral on the whole real line. This approach is equivalent to 
compute R[bj(x)x[-L,L]{\\x\\)](t k ,6 k ), where 

f 1 if — L < r < L 



is the characteristic function of the set [-L, L] for some L > 0. 

Moreover, in general, we can multiply bj for a window function w, where 
w is such that 



Possible choices of w are: 

• the characteristic function of a compact set w(x) = X[-l,l]{\\%\\)', 

• the Gaussian function w(x) = e~ e2 ^ 2 ; 

• the cosine window w(x) = cos 2la[-l,l] ll x ll j 

This approach can be interpreted also as substituting the operator R with 
another operator, say R w , defined by 



Note that, since R is a linear operator, also R w is so. Indeed, for all /, g 
functions, for all a, (5 constants 





otherwise 




R w [f] = R[fw], 



for all / : R 2 



-> R. 



R w [af + Pg] = R[(af + f3g)w] = R[afw + /3gw] = 
= aR[fw] + /3R[gw] = 
= aR w [f]+(3R w [g}. 
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Then if we approximate f k = R[f](t k ,9 k ) with f Wtk = R w [f](t k ,9 k ), we can 
consider the interpolation problem 

fk - fw,k = R w [s](t k ,9 k ) Vfc = l,...,n. 
By linearity of R w 

n 

fk ~ YjCjR^i^k^k) V/C = 1, . . . ,11, 

3=1 

that leads us to the linear system A w c = f, where (A w ) k j = R w [bj](t k , 6 k ) = 
R[bjw\(t k ,9 k ). 

We notice that for all k, the difference between f k and f W}k is bounded by 



\fw,k — fk\ — \R w [f]{tk, Ok) — R[f]{t k ,6 k )\ — 

= \f f(x(s))w(x(s))ds- f f(x(s))ds\< f \f\\w 

\JR JR I JR 

< \\W - l||oo||/||ii(K2) 



lids < 



Thus, for w ->• 1, \f Wik - f k \ ->• but also R[bjw](t k ,9 k ) ->■ 
and this quantity can be infinity. For to — >■ 0, 9 k ) — > but the 

difference 



\fw,k ~ fk\^ J f(x(s))ds 



< ii ; i| L i. 



Before starting on computing y4 w consider the following example. Let K 
be the inverse multiquadric kernel given by 

K(x,y) = . 

v 1 + \\ x - y\\ 2 

As before we can just consider R[k(y)] = R[K(0,y)] because of the relation 
k(y — x) = K(x, y) and the shift property of the Radon transform. What we 
obtain is 

fl[%)iM)=X Vl+ ; 2+ ^=+oo, 

that means that in this case, not only bj(x) is not integrable on some line (as 
in the Gaussian kernel case), but even K{x,y) is not integrable on any line 
l t fi- In this case we have to consider a further regularization of the integral. 

The remedy we adopt is to multiply the function K itself by a window 
function w such that R y [K(x, y)w](t, 9) exists finite for all (t, 6). In choosing 
the function w, we consider that we would like to still use the shift property 
of the Radon transform, therefore we take w of the form w = w(x,y) = 
w(\\x — y\\) so that, if we set now k(y) = K(0,y)w(0,y), it is still true that 
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k(y — x) = K(x, y)w(x, y). Moreover we will choose w to be positive definite, 
in this way, since the product of positive definite function is positive definite, 
also K ■ w is so. 

Notice that this kind of regularization does not correspond, as in the 
first case, to replace the operator R with another operator. In fact now the 
function w depends on x. What we are doing now is simply substituting the 
positive definite kernel K(x, y) with another positive definite kernel given by 
K(x,y)w(x,y) that is integrable on every line in the plane (2/1 , 2/2) - 

5.1.2 Regularization by truncation 

We first consider the regularization of the Gaussian reconstruction problem 
using w L (x) = X[-l,l](\\x\\) as window function. With this choice R[bj](tk, Ok) 
is replaced by R L [bj](t k , 9k) = R[bjW L ](t k , k ), where bj(x) = R y [K(x, y)](tj, 0j) 
and K the Gaussian kernel. 

In applications is useful to use kernels depending on a shape parameter e 
so that choosing suitable values of e > one can obtain system matrix with 
a better condition number. We will consider 



K(x,y) = e 



Basis bj then becomes 



Indeed 



bj(x) = ^V^fe-W 2 . 



£ 



(5.2) 



-e 
e 



R[k(y)}= ( e- s2( - t2+s2 Us = e-^ 2 / e~ s2s2 ds 
Jr Jr 

Components of matrix A L are given by akj = RL[bj]{tk,0k), so we have 

R L [b j ](r lV )= [ ^e-^-^Xi-LMMmds, 

where x(s) = (r cos ip — s sin tp, r sin ip + s cos p) and v = (cos p, sin p). If we 
set again 

a = sin (p — 6) b = t — r cos (p — 6) (5.3) 

we have 



= ^[^e-^ 2 ds, 

£ J-VL 2 -r 2 

where we are assuming L 3> so that \r\ < L. Now we distinguish two cases: 
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1 . if a = then 

R L [bj}( r , <p) = ^ e^ h2 ds = v^ e -^2v^H; 

2. if a 7^ we set u = e(as + 6) so that 

e z a J Cl 

where C\ = e(—^/L 2 — r 2 +b) and C2 = e(\/L 2 — r 2 +b) and J^ 2 e - " 2 cfot = 
^erf(c2) — erf(ci), where erf is the usual error function 

2 f X 2 

erf(x) = — p= / e~" du. 
a/7t Jo 

In conclusion we have to solve the linear system A L c = f with components 
of A L given by 




where 

a = sin (6k — 9j) b — tj — tk cos (6k — 6j) 

Cl = £ (-^/L 2 -t 2 k + b) c 2 = e(^L 2 -t 2 + b). 

Solving this system we obtain c and we can then evaluate the solution s by 

n 

S(x) =Y,Cjbj(x) 

i=i 

with bj given by (5.2). Figure 5.1 shows the results of applying this method 
to the crescent-shaped phantom for suitable values of e and L in the case 
of parallel beam geometry, where samples of the Radon transform of / are 
taken at angles 6 P = p-rr/N, i = 1, . . . , N — 1 and t q = q/M, q = —M, . . . , M. 
We also observe that without using the shape parameter e, i.e. using e — 1, 
the matrix A L would have been highly ill-conditioned and the result very 
different. 
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(c) (d) 

Figure 5.1: Reconstruction of the crescent-shaped phantom with Gaussian 
kernel and truncation regularization : (a) Original phantom; (b) e = 1, 
L = 10, K\A) = 4.04307- 10" 28 ; (c) e = 25, L = 10, K\A) = 2.3734- 10" 5 ; 
(d) e = 50, L = 10, ki\A) = 4.4333 • 10" 5 . 
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5.1.3 Regularization by Gaussian filtering 

Consider again the Gaussian kernel K(x,y) = e~ e2 W x ~y\\ 2 and the associated 
basis of the space S, defined in (5.2). In this section we will use another 
window function in order to regularize the integral R[bj], i.e. we will multiply 
the function bj by another Gaussian function 

Wv {x) = e-" a ll*ll a . 

In other words we will consider the operator R v given by R„[f] = R[fw v ] 
instead of the classical Radon transform. We have 



R v [b j ](r,<p) = [ bJx(s))e-» 2 ^ 2 ds = [ e -eHas+bf & -u^W) ds 
JR £ Jr 

where a, b are defined by (5.3). Then, 



exp (-[(eV + u 2 )s 2 + 2abse 2 ]) ds 



a/7T / 2 2 2 2 a 2 b 2 £ 4 

= — exp -e 2 b 2 - v 2 r 2 + 



exp. 



/ 



exp 



TT 



Ve 2 a 2 + v 2 s + 



a 2 e 2 + v 2 
abe 2 



exp 



-v 2 [ v 2 + 



\Je 2 a 2 + v 2 

e 2 b 2 
a 2 e 2 + v A 



ds 



£\Ja 2 e 2 + v 2 

We have now two options: 

1. The regularization R u \bj](t k ,9 k ) for all values of k,j, that leads to a 
linear system with matrix A\ whose components are 



a k,j — 



tt exp 



2 ( 2 , e 2 b 2 



a 2 e 2 + v 2 



2. The regularization Rv[bj](tk,0k) only for those values of k,j for which 
bj has not finite Radon transform, i.e. only when a = 0, while for 
a^Owe consider the usual Radon transform. This corresponds to the 



matrix A v 2 whose elements are 



ahj - < 



7rexp \-{v 2 r 2 + e 2 b 2 )] 



if a ^ 



if a = 0. 



ev 
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Numerical experiments (Figures 5.2 and 5.3) show that the first option gives 
better results, provided that the value of e is relatively big (« 30) and value 
of v is quite small (« 0.5) so that the condition number of the matrix A\ is 
small. 




(a) Option 1 (b) Option 2 

Figure 5.2: Gaussian reconstruction, parallel beam geometry data 




(a) Option 1 (b) Option 2 

Figure 5.3: Gaussian reconstruction, scattered data 
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5.2 Inverse multiquadrics reconstruction 

We now consider the same reconstruction problem by using the inverse mul- 
tiquadrics as kernel function: 

K(x,y) 



\Jl + e 2 \\x — y\\ 2 

As we saw in section 5.1.1 K does not admit finite Radon transform and so we 
have to multiply K by another window function w of the form w = w(\\x— y\\), 
so that R y [K(x, y)w](t, 9) < oo for all t, 9. 

The window function we consider is the characteristic function of a com- 
pact set: 

w = X[-L,L](||z-y||), £»o. 
As in the Gaussian case we use the shift property of the Radon transform 
and we first compute 

R L [k(y)](t,e) = R[K(0,y) X[ _ LM (\\y\\)}(t,9) = 



= L 2 ,- X[-L,L](VF + ^) ( fa 

JMyfl + e 2 (t 2 + s 2 ) 

rVL 2 -t 2 I 

= / . — ds 

J — - 



fVi 2 -* 2 l 
-VTJ^ 1 y/c 2 + (es) 2 

where c = VT+ e 2 t 2 and we assume \t\ < L. We then apply the substitution 
u = es in the integral 

1 r cVL^ i , \2 . , /u\l e ^^ 2 . / lL 2 -t 2 \ 



£ 



du = -asinh 



where 



= e asinh ^Vi + ^ 2 / 

(5.4) 



asinh(x) = log (x + VT+ x 2 ), 
so we can also write integral (5.4) as 

e^JJ^I 2 ^eVL^fi 



- [asinh (") 



- [log(u + Vc 2 + "« 2 )] V 
= 2 - (tog {eVL 2 - t 2 + Vl + e 2 L 2 ) - * log (1 + e 2 * 2 )) . 



We conclude that the basis associated to the kernel Kw is 
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where, as usual, Vj = (cos 9 j, sin 9 j). 

We can now compute the matrix A = (cik,j)k,j=v 



f 2 / L 2 — (as + b) 2 \ 
a kJ = R[bj(x)](r, <p) = ^ -asinh ( eJ Y+e^as + b) 2 ) x ^ L > L ^ as + &) ds ' 

where a, b are defined as usual and r = t^, if = 9 k . We observe that if a = 0, 

then 

2 . / \1? 



a k!j = -asinh ^ _j ^ X[ _ L , L] (6) <fe = +oo 

when |6| < L, that is our case. So we have to consider a further regularization 
of i?. We choose to truncate, i.e. we compute 

OfcJ = ifo[&,-(aO](r,¥>) = i%(x)x[_H,ff](||z||)], H > 0. 

As for the Gaussian case, we have two options: consider the regularization 
Rh for all values of k and j or use it only when 9k = 0j i.e. when a = 0. 
As before we consider the hrst option since the resulting matrix has a better 
condition number. Thus, we obtain 

a k ,j ^//^( ^ + " £ ^^ 2 2 ) x [ -L M (as + b)x [ -H,H ] (^T^)ds 
that in the case a = becomes 



4 / L 2 — b 2 \ 
a kJ = -asinh J ^H 2 - r 2 , 

provided that \b\ < L and \r\ < H. In the case o ^ we consider the 
substitution u = e(as + b) that leads us to the integral 

(5.5) 




Ci = emax (-L, —\a\\/H 2 — r 2 + b) C2 = emin (L, \a\ v 7 H 2 — r 2 + b). 

All what we have to do now is to compute integrals (5.5). The computation 
of these integrals can be found in appendix A. 

Applying this method again to the crescent-shaped phantom, choosing 
e = 30, H = L = 20max|tj| we obtain the reconstruction shown in Figure 
5.4. We observe that we have acceptable reconstruction only with paral- 
lel beam geometry data, to obtain good results also with scattered data we 
should consider another window function instead of the characteristic func- 
tion. 
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(a) Parallel beam geometry 



(b) Scattered data 



Figure 5.4: Inverse multiquadric reconstruction of the crescent-shaped phan- 
tom using e = 30, H = L = 20. 

5.3 Multiquadrics reconstruction 

We now consider the multiquadric kernel 



As in the case of inverse multiquadrics, K is not integrable on any line in 
the (2/1,2/2) plane. The approach we follow to regularize the problem is to 
consider a Gaussian weighting function for computing both basis functions 
bj and the matrix A. 

We start with the Gaussian-filtered basis 



We recall that this operation corresponds to use kernel K(x, y) = K(x, y)e~ £2 ^ x ~ y ^ 
in place of K. Proceeding as usual, we first compute bj(x) then the coeffi- 
cients a,j t k- 



K(x, y) = ^Jl + p 2 \\x-y\\ 2 , p>0. 



bj(x) = R y {K(x,y)e- £2llx - yll2 }(t j ,e j ), e > 0. 



Now, 




p 2 as. 
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Setting c = \J~jp +t 2 and then integrating by parts, we obtain 
1 



bj(0) = 2pe^ \c 2 

f+OO 



+oo 



/ 

JO 



4 sinh (2*) + g + 
^ sinh (2s) + |) e - £2 ( c2 +- 2 ) (-2s 2 s)e-^ c2+s ^ ds 



c 2 pe p' 2 



^sinh(2s) + s)e- £2 ( c2+s2 ) 



+0O 



+ 



e 2 „ +0 o / I \ 

+ 2c 2 pe» 2 £ 2 J ( -sinh (2s) + sj 

e 2 e 2 

= c 2 pe^h + 2c 2 pe^e 2 I 2 , 

where 

ii = lim 



s=0 



se~^ + ^ ds = 



r— »+oo 



sinh(2r) + r) e^ 2 ^ 2 ) 



e- £2c2 (^sinh(0) + 0) = 



and 



1 r+°° o,o 9n r+oo 

h = - / s sinh (2s)e" e (c +s > ds + / 

2 Jo Jo 



sV^ c2+ * 2 > ds 



2 >— 1 /— -£ 2 C 2 

7T 1 \/ 7T e .-2 

e . 



1 _ e 2 c 2V^ e ° . - £ 2 e 

_p — L p fc L . . 

2 2 e 3 4 £ 3 4 e 3 



We conclude that 



bj(0) 



\pnp ( 1 



2e \p 



2 \ „-e 2 t 2 +e- 



and so 



6 J (x) = ^(l + (t J -x.t, i ) 2 ) 



-e 2 (tj-a;-ii ; ,') 2 +e 2 



= (cos#j, sin^j). 



In order to compute the matrix A, we consider R[bj(x)](r, <p) . Setting 
x s = (r cos 6 — s sin 0, r sin 6> + s cos 6>) , we have 

a kj = RMx)](r, tp ) = ^^(\ + (t-x 8 - v) 2 ^j e -e 2 it- Xs -v ?+ e- 2 ds = 



2e f 



/ 

■m 



p e 



2 p -e 2 (as+bf dg + 



f (as + b) 
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where a = sin (tp — 9) and b = t — r cos (ip — 6). If a = this integral is 
infinite. To avoid this case we consider the regularization R v of R, that is 

RAbj(x)](r, <p) = [ f 1 + (t - x s ■ v)A e -*t-*.«r e -*M> ds = 

ZE Jr \p ) 



" 2 / ( L + (GS + b) A e -^M)M^) ds 

JR\P 2 ) 



2 e-^ 2 ds + 



I {as + b) 

Jr 



2 e^ cs+d y 2 ds 



where 



c = Vv 2 + sW, d=—, C £v = ^exp(e- 2 + d 2 -e 2 b 2 -v 2 r 2 ). 

c 2e 



Since 



/ (as + b) 2 e-( cs+d f ds = ^(a 2 (2d 2 + I) - Aabcd + 26 2 c 2 ), 
Jr 2 c r 



we conclude that 



vrexp(^- ^^-//V 2 ) M pa 2 (^ 2 + £ 2 a 2 ) + 26V 



2eV^ 2 + e 2 a 2 



+ 



p 2 (v 2 + e 2 a 2 ) 2 



Figure 5.5 shows the result of using the multiquadric kernel with Gaussian 
filtering for the reconstruction of the crescent-shaped phantom. 



5.4 Compactly supported radial basis func- 
tions 

Another important class of positive definite functions we consider are the 
compactly supported radial basis functions. If a function as compact support, 
then it is automatically strictly positive definite and it is strictly positive 
definite on M d only for a fixed maximum value of the dimension d. Moreover 
one can show that there not exist compactly supported radial functions that 
are strictly conditionally positive definite of order m > (see [23] for more 
informations about compactly supported radial functions). 



Wendland's compactly supported functions 

A popular family of compactly supported functions was introduced by Wend- 
land [22]. Wendland starts with the truncated power function (pi(r) = 
(1 — r)' + , which is striclty positive definite and radial on R d for d < 21 — 1, 
and then applies repeatedly the integral operator X defined as follow: 
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(a) Parallel beam geometry 



(b) Scattered data 



Figure 5.5: Multiquadric reconstruction of the crescent-shaped phantom us- 
ing = 1, e = 30, v = 0.8. 

Definition 13. Let ip : [0, oo) — > M. such that t<p(t) 6 L^O, oo), then we 
dehne 



lip(r) = / t(p(t) dt, r > 0. 

J r 

We can now define the Wendland's compactly supported functions: 
Definition 14. With <pi(r) = (1 — r) l + , we define 



Example 9. The explicit representation of ip<i t k for k — 0, 1, 2, 3 are: 



<Pd,o(r) = (1 -r) l + , 

ip d>1 (r) = (l-r) l + [(l + l)r + i\, 

^, 2 (r) = (1 - r)' + [(/ 2 + 4/ + 3)r 2 + (3/ + 6)r + 3], 

^ 3 (r) = (1 - r)^[(/ 3 + 9l 2 + 23/ + 15)r 3 + (6/ 2 + 36/ + 45)r 2 + (15/ + 45)r + 15], 



where / = [d/2\ + k + 1 and equalities are up to a multiplicative constant. 

We observe that all the functions in Example 9 are compactly supported 
and have a polynomial representation on their support. This is true in general 
as stated in the following 

Theorem 5.4.1. The functions ipd,k are strictly positive definite and radial 
on M. d and are of the form 



'OO 
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where pd,k is a polynomial of degree [d/2\ +3fc + l. Moreover ipd,k £ C 2k (R d ) 
are unique up to a constant factor and the polynomial degree is minimal for 
given space dimension d and smoothness 2k. 

The proof of this theorem can be found in [22] . 




Figure 5.6: Plot of Wendland's compactly supported functions. 



5.4.1 Compactly supported kernel reconstruction 

Compactly supported radial basis functions can be used as well in image re- 
construction: if <p(r) is a compactly supported function, one sets K(x,y) = 
<p(s\\x — y\\), e > 0, and uses K as kernel function. The property of compact 
support can be useful in the computation of the Radon transform, in partic- 
ular if ip is compactly supported and of class at least C° on its support, then 
the Radon transform Rf is well defined, indeed 



\Rf\ = 


/ f(x 8 )ds 




/ f(x s ) ds 


< 1 






JSupp(f) 


JSupp(f) 



ds < 



< sup \f{x)\ ■ m(Supp(f)) < oo. 
xeSupp(f) 



For example Wendland's and Wu's compactly supported functions, having 
a polynomial representation on their domain, admit finite Radon transform. 
This means that for this class of functions the basis bj(x) = R y [K(x, y)](tj, 9j) 
is well defined. 
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However the fact that the kernel K is compactly supported does not imply 
that bj(x) has finite Radon transform and then matrix A = (a^) can have 
non finite entries as shown in the following example. 

Example fO. Consider the Wendland's function (f2,i(r) = (1 — r) 4 (4r + f) 
and set 

K{x,y) = ^2,\{£\\x - y\\) = (1 - e\\x - y||) 4 (4e||x - y\\ + f) 

The support {(x,y) : e\\x — y\\ < 1} of K is compact. 

We compute bj(x) using the shift property of the Radon transform: 

R y [K(0,y)](t,9) = [ (f - eVt 2 + s 2 )1_(4ev / t 2 + s 2 + f ) ds = 
= [ 4e(l - eVt 2 + s^y/PTs* ds+ 

+ / (f-eVt 2 + s 2 ) 4 ds, 



thus 



R»[K(0,y)](t,9) 



h + h if \t\ < 



if \t\ > 



with 



h= I . 4eVt 2 + s 2 (l - eVt 2 + s 2 ) 4 ds = 

= - ipi(et)acosh (—rn J — p2(st)Vl — £ 2 t 2 



where 



Viir) = -(5r 4 + 36r 2 + 8), 



Mr) = ^(437r 4 + 360r 2 -8) 



acosh(r) = log(r + Vr 2 — f). 
While the second integral is 

h= I , (1 - eV^T^ 2 ) 4 ds = 



-p 3 (rf)acosh (— !— J + p 4 (et)v / l — e 2 t 2 
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1 

15 



(16r 4 + 83r 2 + 6). 



with 

Ps(r) = r 2 (3r 2 + 4), p 4 (r) 
Then we obtain 

h + h = - £ (pi - p 3 )(^)acosh ^-j^ + (p 4 - p 2 )(et)v / l-e 2 t 2 



V(e 2 t 2 + 6)acosh f ^] - 1 (81e 4 t 4 + 28e 2 t 2 - AWl-eH 2 
\e\t\ / 6 



We observe that I\ and 7 2 are not well defined for t = 0, in this partic- 
ular case it is easy to prove that the value of the integral is ^, moreover 
lim^o (h + I2) = y so we can define the continuous function 



9e(t) 

then bj(x) is given by 
bj(x) = 



' h + h if e\t\ < 1, t ^ 
2 

- 3e 



if t = 0, 



5 £ (ij — a; • Wj) if eft,- — x • t> j\ < 1 







if e\tj — x ■ Vj\ > 1. 



For fixed x, it^fi^x, 0) as function of (t, 0) has compact support, but 
if we consider fixed values (tj,9j), then R y [K(x, y)](t, 9) = bj(x) as function 
of x has {x e M 2 : — x ■ vj\ < 1} as support, that is the strip between 
lines e(tj — x ■ Vj) = ±1 that is not limited and thus not compact. 

If we now compute cikj = R[bj(x)](tk,0k) we see that this quantity can 
be infinity: 

a k j = R[b(x)](r, <p) = / b(x s )ds, 

where x s = (r cos ip— ssinip, r sinip+s cos ip). We set as usual t—x s -v = as+b, 
then 



a k,j — / 



'e\as+b\<l 

If a 7^ we can set u = s(as + b) obtaining 



g £ (as + b) ds. 



• - k L + 6,acosh (r) - s (81t * 4 + - 2 - 4) ^ 



du 



9 7T 

112 e^a 
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But if a = we have 



a k,j — 




Je\b\<l 



ds 



that is if s\b\ > 1 but is oo if e\b\ < 1. 

In order to have a matrix A with all finite entries, we again consider 
the regularization R w of R for some weighting function w. Working with 
compactly supported radial basis functions it's natural to use another com- 
pactly supported function as weighting function, in this way we are sure that 
R w [bj] = R[bjw] is finite (indeed bw is compactly supported) and it is possible 
to compute analytically the value of a^j. 

We consider the following case: 



K(x, y) = ip 2fi = (1 - e\\x - y||)+, 
w(x) = (1 - v 2 \\x\\ 2 ) + , v > 0. 



e > 0, 



We start by computing 



R[K(0,y)](t,e)= [ {l-eVWT^) 2 + ds= [ 

_ ' g(t) if 1*1 < \ 



(1 - sVt 2 + s 2 ) 2 ds 



if \t\ > - 



e 



where 




= < 



' 2 r yi -sH 2 

£ [ 3 

2 

. 3^ 




if t = 



We conclude that 
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As in Example 10, it is possible to show that akj = oo if a = and 
e\b\ < 1, where a = sin (0*. — #j) and b — tj — tk cos (0*, — fy), so we introduce 
the regularization of i? by multiplication for the weighting function w = 
(1 — z/ 2 ||x|| 2 ) + and we have 

ak,j = Rwibjix^tk^k). 

Using simpler notation, 

R[wb(x)](r,(p)= f b(x s )(l-u 2 \\x s \\ 2 ) + ds. 

JR 

The computation of this integral can be found in appendix B. 



5.5 Scaled problem 

The Hermite-Birkhoff reconstruction problem can be expressed as finding a 
function seS satisfying 

\ j (s) = \ j (f) Vj = l,...,n, (5.6) 

where Ai, . . . , A„ are linearly independent linear operators and 

r „ \ 

S = <Y,CjXjK(-,y) ■ Cj e R, ifposit ive definite kernel > . 
b=i J 
By linearity equation 5.6 can be written as a linear system Ac = /|a, where 
the elements of the matrix A are given by a^j = \%\V[K(x, y)] and (/|a)j = 

W). 

Since the matrix A can be highly ill-conditioned, it can be convenient to 
consider the scaled problem ([10], [11]). For h > the scaled reconstruction 
problem is 

\ J (s h (h.)) = X J (f(h.)) Vj = l,...,n, 

where 

s h eS h = jfj CjA^(-, %) : \ Cj e mJ . 

In this way one obtains the linear system A h c — f h \\ given by 

4j = \l\ V j[K(hx,hy)], (/ fc U)i = ^j(f(h-)). 

In the case of image reconstruction the operator Xj represents the Radon 
transform evaluated at point (tj,9j). Thus, in order to compute a\- and /j 1 
one has to understand which is the relationship between the Radon transform 
of a function / and the Radon transform of the scaled function f(h-). This 
relationship is given by the following 
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Theorem 5.5.1 (Dilatation-property of the Radon transform). Let f : M 2 — > 
R be such that R[f(x)](t,6) = g(t,9), then for all h > 

R[f(hx)](t,e) = p(ht,e). 

Proof. 

R[f(hx)](t,6) = f f (ht cos 9 - hs sin 9, ht sin 9 + hs cos 9) ds = 

f dr 
= / f (ht cos9 — rsin9,htsm9 + r cos 9) — = 
Jr h 

= l l R[f{x)](th,9). 

□ 

Thanks to this property it's possible to compute and bj(x) = 

\?K(x,hy): 

ti = \l{f{hx))= l -R[f{x)]{ht kl 9 k ) 
b h ] {x)= l J Ry[K(. 1 y))(ht ]1 9 3 ) 

a h kJ = \%\M[K(hx,hy)] = \\ x k [b){hx)} = ^R[bj(x)]{ht k ,0 k ) 

We note that in order to compute f k we need to know the Radon trans- 
form of the unknown function / at points (ht k , 9 k ) that is not possible if the 
X-ray machine gives us data only at points (t k ,9k). However, for our aim, 
we assume to know the analytical expression of R[f]. 

The solution of the scaled reconstruction problem is 

n 

s h (hx) = Y,Cjtf{hx), 

3=1 

with c = (ci, . . . , c n ) T solution of the linear system A h c = f h . It is important 
to notice that thanks to the dilatation-property 5.5.1 we don't have to com- 
pute the Radon transform for every different value of h, but we only need 
to compute it in the case h — 1, multiply for ^ and then scaling evaluation 
points (t k ,0 k ) to (ht k ,O k ). 



Chapter 6 
Numerical results 



In the previous chapters we saw some theoretical tools that can be used to 
obtain the value of a function / : M 2 — > M starting from a sampling of its 
Radon transform. 

In chapter 2 we studied the continuous problem and we found an analyt- 
ical inversion formula for the Radon transform: the back projection formula. 
Then, in order to use this formula in real applications, we introduced the 
process of linear filtering and interpolation. 

In chapter 3 and 4 we followed a different approach: starting form the 
discrete problem for finding an approximation of a function, belonging to a 
particular finite dimensional space of functions, such that its Radon trans- 
form coincides with the measured Radon transform of the unknown function 
/. The Kaczmarz's method consider pixel basis functions to determine an 
approximation of /, while kernel-based methods use positive definite func- 
tions to generate a functions space where to find a solution. We also saw 
that in this second case the problem need some kind of regularization so that 
the Radon transform of the kernel functions is well defined. 

What we want to do now is to compare all these methods from a numerical 
point of view, studying the behavior of the solution and the approximation 
error in function of the data and the parameters involved in the algorithms. 

6.1 Optimal parameters 

In chapter 5 we introduced a regularization technique for solving the Hermite- 
Birkhoff interpolation problem of image reconstruction using kernel based 
methods. In particular, the original problem was to find s = J2]=i c j^jK(-, y) 
such that Xkf = XkS for all k — 1, . . . , n, where X^g = R[g(x)](tk, 6k) and K 
is a positive definite kernel. By linearity this problem is equivalent to solve 
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the linear system 

n 

x kf = ^2c j \ x k \ y j K(x,y), k = l,...,n. 
3=1 

The main problem found in applying this method is that the Radon trans- 
form X k X^K(x, y) or X V jK(x, y) can be infinity. 

The solution we adopted was to consider kernel functions K such that 
bj(x) = X y jK(x,y) is well defined (e.g. multiplying any kernel K(x,y) for 
a suitable function 4>(\\x — y\\)) and to substitute operator R with another 
linear operator so that, when computing matrix A = (a^-) = (Xkbj(x)), we 
have a k j < oo for all k, j — 1, . . . , n. 

In our discussion we chose operator R w defined by R w [g] = R[gw] where 
w is an appropriate window function. 

Both kernel K and window function w depend on one or more parameters; 
in this section we will discuss, thanks to numerical experiment, how these 
parameters influence the quality of the reconstructed image. In order to do 
that we will apply our methods on predefined phantoms and by varying a 
parameter we will see the behavior of the solution. The error measure we 
will use to determine the quality of the result is the root mean square error 

RMSE = \ ^~ lV — , 

V m 

where m is the dimension of the image, Xi and Xi the gray scale value of pixel 

1 of the original and reconstructed image respectively. More the RMSE is 
close to zero, more the solution will be considered accurate. 

6.1.1 Window function parameters 

We begin our analysis considering the parameter that influence operator R w 
and the window functions introduced in sections 5.2 and 5.1, i.e. w v (x) = 
exp(-// 2 ||z|| 2 ) and w L = X[-l,l](\\ x \\)- 

Let us start with the truncated inverse multiquadric kernel 

K(x, y) = (l + e 2 \\x - y|| 2 )- 1/2 X[- " v\\) 

and the characteristic window function Wl 2 = X[-l 2 ,l 2 ](.\\ x \\)i w hh L 1 ,L 2 > 

2 max |tj| (as we saw in section 5.2). 

Varying the parameter L 2 for fixed values of e, Li and data 1 {(tj, 6 , J )}™ =1 
and applying this reconstruction technique to three different phantoms we 

Mn this chapter we will always consider the parallel beam geometry as acquisition 
method of data. 
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can see that the RMSE presents a minimum for a particular value L 2opt 
(Figure 6.1). This optimal value is influenced by e, Li, n (in particular if 
the number of data increases, also L 2opt increases - cfr. Figures 6.2(a) and 
6.2(b)). We also notice that the RMSE decrease very rapidly for L 2 < L 2opt 
but for L 2 > L 2opt the RMSE is increasing with a very small rate, so one 
should choose L 2 in a way to be sure that L 2 > L 2opt . 




(a) Crescent-shaped phantom (b) Bull's eye phantom 

Figure 6.1: RMSE of inverse multiquadric reconstruction in function of 
parameter L 2 , with fixed number of samples N = 30, M = 20, K = 256 and 
parameters e = 30, L\ = 50. 

Another advantage in choosing L 2 large is that the condition number of 
the matrix A is smaller (see Figure 6.3). 

However the value of L 2 does not determine so drastically the behavior 
of the solution, whose quality remains acceptable for large values of L 2 . 

More interesting is the case of the Gaussian window function w(x) = 

2 II II 2 

e ~ v IfII _ Also in this case there exists an optimal value f opi such that RMSE 
is minimum. But now for v > v opt the RMSE increases with a fast rate and 
so the quality of the reconstruction becomes worse (see for example Figures 
6.4(c) and 6.4(d)). 

The value v opt depends on the phantom used, i.e. on data. This is not 
surprising because we know that the approximation error depends on \f Wt k — 
fk\ (see section 5.1.1). On the other hand v opt has only small variation w.r.t. 
the changing of other shape parameters (e.g. is independent on e in the case 
of Gaussian kernel - see Figures 6.4(a) and 6.4(b)). 

The fact that the RMSE is increasing for v > u opt can be explained 
considering the condition number k(A) of the system matrix A. In fact, 
k(A) increases with v (see Figure 6.5 where the reciprocal of k(A) is plotted 
in function of v). The quantity k(A) = k\(A) is the 1-norm condition number 
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Figure 6.2: RMSE of inverse multiquadric reconstruction as a function of 
the parameter L 2 , with fixed output dimension K — 64 and parameters 
e = 30, L\ = 10 for the crescent-shaped phantom. 
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(a) Parameters and data as in Figure 6.2(a) (b) Parameters and data as in Figure 6.3(b) 



Figure 6.3: Reciprocal of the condition number of A for inverse multiquadric 
reconstruction as a function of the parameter L 2 
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(a) 



(b) 





(c) 



(d) 



Figure 6.4: RMSE of Gaussian and multiquadric reconstruction as a func- 
tion of the parameter v. (a) Gaussian kernel, bull's eye phantom, N = 
50, M = 40, K = 64, e = 30; (b) Gaussian kernel, bull's eye phantom, 
N = 50, M = 40, K = 64, e = 60; (c) Gaussian kernel, Shepp-Logan 
phantom, iV = 30, M = 20, K = 256, £ = 50 (d) Multiquadric kernel, 
crescent-shaped phantom, TV = 30, M = 20, K = 64, p = 1, e = 30. 
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of the matrix A. Its inverse is estimated using the MATLAB function rcond 
(see [8]). 
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Figure 6.5: Reciprocal of the condition number of A for Gaussian kernel 
reconstruction as a function of the parameter v. (a) Parameters and data as 
in Figure 6.4(a); (b) Parameters and data as in Figure 6.4(c). 



The case of multiquadrics K(x,y) = a/1 + p 2 ||a; — 2/ 1 1 e y " with Gaus- 
sian window function is similar to the Gaussian kernel case, provided that p 
is small enough, as explained in the next paragraph. 

At last we consider the case of compactly supported kernel. Let K = 
(1 — e\\x — y\\)\ and w(x) = (1 — // 2 ||z|| 2 ) + . Assume e w 1, then it turns out 
that the RMSE is minimal for v (Figure 6.6). For small values of u also 
the condition number of A is larger (Figure 6.7), so it is convenient to use 
v « 0. 

We observe that using small values of v corresponds to use a compactly 
supported window function w with wide support, in this way one loses less 
information when filters basis bj with w. 

We can use the information that the approximation in optimal for v w 
to simplify the expression of the matrix A. Indeed, when a ^ 0, if v — > 0, 
then a,k,j — > and we can use this simpler expression of atj instead of (.3) 
(see appendix B). This option is equivalent to consider the regularization R w 
only when a = 0, while using the original operator R when a ^ (that is 
what we called option 2 in the section 5.1.3). Using this second option the 
behavior of RMSE becomes more regular (see Figure 6.8). 

Finally we observe that using more data, one should use a bigger value 
of v, as shown in Figure 6.9. 
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(a) Crescent-shaped phantom 



(b) Bull's eye phantom 



Figure 6.6: RMSE of compactly supported reconstruction as a function of 
the parameter v. N = 30, M = 20, K = 64, e = 1.1. 
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(a) Parameters and data as in Figure 6.6 (b) Parameters and data as in Figure 6.6 



Figure 6.7: Reciprocal of the condition number of A for compactly supported 
kernel reconstruction as a function of the parameter v. 
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(a) Crescent-shaped phantom (b) Bull's eye phantom 

Figure 6.8: RMSE of compactly supported reconstruction as a function of 
the parameter v. Regularization only for a = 0. Parameters and data as in 
Figure 6.6. 
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Figure 6.9: RMSE of compactly supported reconstruction as a function of 
the parameter v for the Bull's eye phantom with N = 50, M = 40, K = 
64, e = 1.1. 
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6.1.2 Kernel shape parameters 

We now consider the second main shape parameter involved in kernel meth- 
ods, i.e. the kernel shape parameter e. We will assume to work with an opti- 
mal window function parameter (as discussed in the previous paragraph) . We 
observe that in the following cases the behavior of the RMSE as a function 
of e is similar to that of the RMSE as a function of v. 

• Inverse multiquadrics. There is an optimal value e opt such that RMSE 
is minimum, but for e > e opt the RMSE increase slowly (see Figures 
6.10(c) and 6.10(d)); 

• Gaussian. There is an optimal value e opt such that RMSE is minimum 
and for e > e opt the RMSE increases fast. Moreover e opt increases if 
the the total number of data increases (see Figure 6.10(a) and 6.10(b)). 

The difference with the window parameter is that now the condition num- 
ber k(A) is smaller for big value of e (equivalently the reciprocal fc(^4) _1 
increases with e - Figure 6.11). 

Another parameter involved in the inverse multiquadric reconstruction 
is L\. We observe that under a certain threshold the behavior of the error 
as a function of L\ is chaotic, but over this threshold the RMSE remains 
almost constant. The reciprocal of the condition number is instead decreasing 
(Figure 6.12). 

Consider now the Gaussian-multiquadric kernel 

K(x,y) = ^l + p 2 \\x-y\\ 2 e- £ ^ x - y ^ 

. In this case we have two different shape parameters p and e. As in the 
Gaussian case, large values of e generate matrix with a better condition 
number (p fixed). For fixed e » /) we have K(x,y) ~ e~ £ \\ x ~ y W a kernel 
similar to the Gaussian case, but values p « give a bigger conditioning. 

If p is of the same order of e then k(A) is small but MSE becomes larger. 
So we have again a situation with an optimal value for p and e that arises 
from the trade-off to have a well conditioned matrix and a good approxima- 
tion of the non-regularized reconstruction problem (trade-off principle [20]). 
Moreover, as well as e also the optimal value of p depends on the num- 
ber of data and on the phantom. For example, using the crescent-shaped 
phantom, the value of p opt varies from ps 3 (for N — 20, M — 15) to « 7 
(for N = 50, M — 40); while considering the Shepp-Logan filter we obtain 
Popt < 0.5 (for N = 20, M — 15) and p opt < 2 (for TV = 50, M — 40) (Figure 
6.14). 



86 



CHAPTER 6. NUMERICAL RESULTS 




(a) 



(b) 



0.22 



0.215 



0.21 



0.205 



0.418r 




100 



(c) 



(d) 



Figure 6.10: RMSE of Gaussian and inverse multiquadric reconstruction as 
a function of the parameter e. (a) Gaussian kernel, crescent-shaped phantom, 
N = 30, M = 20, K — 64, v = 0.5; (b) Gaussian kernel, bull's eye phantom, 
N = 50, M = 40, K = 64, v = 0.7; (c) Inverse multiquadric kernel, 
crescent-shaped phantom, TV = 30, M = 20, K = 64, L x = 10, L 2 = 10; (d) 
Inverse multiquadric kernel, Shepp-Logan phantom, TV = 50, M = 40, K — 
64, L x = 10, L 2 = 10; 
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Figure 6.11: Reciprocal of the condition number of A for Gaussian and in- 
verse multiquadric kernel reconstruction as a function of the parameter e. 
(a) Parameters and data as in Figure 6.10(a); (b)Parameters and data as in 
Figure 6.10(d) 




(a) RMSE (b) k(A)- 1 

Figure 6.12: RMSE and k~ 1 (A) for the inverse multiquadric reconstruction 
of the Shepp-Logan phantom as a function of Li with N = 36, M = 25, K — 
60, e = 27, L 2 = 10. 
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(a) RMSE (b) k(A)- 1 



Figure 6.13: RMSE and k 1 (A) for multiquadric reconstruction of the 
crescent-shaped phantom as a function of p with N = 50, M = 40, K — 
64, e = 52, i/ = 0.5. 
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(a) TV = 20, M = 15, e = 20.4 

Figure 6.14: RMSE for multiquadri 
function of p with if = 256, // = 1.3 
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At last we consider again the compactly supported kernel. Varying e we 
see that the RMSE presents a minimum for an optimal value s opt that is 
between 1 and 2 for both the crescent-shaped and the bull's eye phantoms 
(Figures 6.15(a) and 6.15(b)). Considering instead the Shepp-Logan phan- 
tom the RMSE decreases if e increases (Figure 6.15(c)). However, in all 
the cases, the RMSE remains almost constant for e large enough, while the 
reciprocal of the condition number increases with e (Figure 6.16). 




(a) Crescent-shaped phantom (b) Bull's eye phantom (c) Shepp-Logan phantom 



Figure 6.15: RMSE of compactly supported reconstruction as a function of 
the parameter e. N = 30, M = 20, K = 64, v= 10" 6 . 




Figure 6.16: Reciprocal of the condition number of the matrix A of compactly 
supported reconstruction as a function of the parameter e with N = 30, 
M = 20, K = 64, e = 1.1. 



6.1.3 Scale parameter 

Considering the scaled reconstruction problem of section 5.5, one has also to 
consider the behavior of the solution depending on the scale parameter h. 
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Numerical experiments show that there is an optimal value of h such that the 
RMSE is minimum and k~ 1 (A) is maximum (Figure 6.17(a) and 6.17(b)). 
It turns out that, in most of the cases we examined, the optimal value of h 
is h « 1. Thus, in our discussion we will consider h — 1. 




(a) RMSE (b) k-^A) 



Figure 6.17: Gaussian kernel reconstruction of the crescent-shaped phantom 
in function of h. Here e = 50, v = 0.5, N = 30, M = 20, K = 64. 



6.2 Comparison of the methods 

In this section we compare the classical Fourier methods, introduced in chap- 
ter 2, with the kernel-based methods of chapter 5. In this second case we 
assume the use of optimal shape parameters. 

We compare the solutions of different algorithms varying the phantom 
and also testing their behavior when introducing some noise in the data. 
Again we use the RMSE to measure how much the solutions differ from the 
original phantom. 

We start considering the behavior of the RMSE in function of the number 
of the data n (where again data are supposed to be taken using a parallel 
beam geometry). As one would aspect, with both kernel and Fourier-based 
methods, the RMSE decreases when n increases. In particular is interesting 
to notice that in the Fourier reconstruction, the RMSE decreases with an 
exponentially rate and so, for large n, there is no big improvement of the 
solution. For example in the case of the Shepp-Logan phantom, for n > 18090 
there is a variation of the RMSE lower than 5.67 • 10~ 3 (Figure 6.18(a)). 

In the case of kernel reconstruction the use of big amount of data should 
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(c) (d) 

Figure 6.18: Reconstruction of the Shepp-Logan filter using the back projec- 
tion formula with Shepp-Logan filter and linear interpolation, (a) RMSE as 
a function of the number of data (without noise); (b) RMSE as a function 
of the number of data (with Gaussian noise with mean ji = 0.01 and variance 
a = 0.01); (c) Reconstruction without noise with N = 180, M = 200; (d) 
Reconstruction with Gaussian noise with N = 180, M — 200. 
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consider CPU limits 2 . Indeed the matrix A used to compute the coefficient 
c of the solution belongs to the space Mat(n,n), moreover, to evaluate the 
solution on a grid of K x K pixels, we must multiply c for a matrix B e 
Mat(K 2 ,n) representing the the basis functions of the space our solution 
belongs to. Thus, for example, with K = 256, N = 50, M = 40, one obtains 
two (non-sparse) matrix one with (N* (2M + 1)) 2 = 16.402.500 elements and 
the other with K 2 ■ (N * (2M + 1)) = 265.420.800 elements. Furthermore we 
notice that if n increases, also k(A) increases (see Figure 6.20(b)). 

Considering a problem with reasonable values of n and K (e.g. n < 
130, K < 256), we observe that the RMSE of the kernel methods behaves in 
the same way as the Fourier based methods and has a comparable magnitude 
(Figure 6.20(a)). In particular, Figure 6.19 shows how the Gaussian kernel 
reconstruction applied to the crescent-shaped phantom gives better result 
w.r.t. to the back projection formula. 

Thus we conclude that kernel based methods can be useful in the context 
of limited number of available data. 

Comparing the reconstruction of a phantom using different kernels, we see 
that the Gaussian-multiquadric kernel and the Gaussian kernel give similar 
results while the truncated-multiquadric kernel and the compactly supported 
kernel are less accurate (Figure 6.21). 

The biggest problem in the kernel methods, a part from the memory lim- 
itation, is the computational time. Referring to Figure 6.20(c) and 6.20(d), 
we can see how the elapsed time (in sec) during the execution of the algorithm 
grows exponentially with the number of data, while using Fourier techniques 
the time depend linearly on the dimension of the problem. We belive that 
this is due to the implementation in MATLAB of Fourier transform by the 
FFTW algorithm [6]. 

Finally we test our methods after the introduction of noise in the data. 
We hrst notice that for large n the RMSE of the Fourier methods increase 
(Figure 6.18(b)). Considering a kernel based method we observe (in the range 
of acceptable n) a behavior similar to the Fourier case, where the RMSEs 
computed with different methods have the same order of magnitude (Figure 
6.22 and 6.23). 

6.3 Graphical user interface 

In order to test the various methods on different phantoms we developed a 
graphical user interface (GUI) allowing the user to choose the options of the 



2 A11 the computation shown in this chapter are made using a computer with a CPU 
core i5, 2,53 GHz and a RAM of 4 GB. 




Figure 6.19: Comparison of Fourier and Gaussian kernel reconstruction meth- 
ods of the crescent-shaped phantom: (a) Root mean square error; (b) Fourier 
method N = 50, M — 40; (c) Gaussian kernel N = 50, M = 40. 
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Figure 6.20: Comparison of the Fourier and multiquadric kernel reconstruc- 
tion of the Shepp-Logan phantom: (a) Root mean square error; (b) Recip- 
rocal of the condition number of matrix A of the kernel method; (c) Elapsed 
time for Fourier based method; (d) Elapsed time for multiquadric reconstruc- 
tion method. 
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(c) 



Figure 6.21: Comparison of kernel reconstruction methods of the bull's eye 
phantom using different kernel functions: (a) Root mean square error; (b) 
Reconstruction with inverse multiquadrics kernel N = 50, M = 35; (c) 
Reconstruction with Gaussian kernel N = 50, M = 40. 
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Figure 6.22: Comparison of the Fourier and Gaussian kernel reconstruction 
of the crescent-shaped phantom with no noise data and Gaussian noise data 
with 0.001 mean and 0.001 variance: (a) Fourier method; (b) Gaussian kernel 
method. 
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(a) (b) 

Figure 6.23: Comparison of the Fourier and Gaussian kernel reconstruction 
of the crescent-shaped phantom with Gaussian noise data (mean fi = 0.001 
and variance a = 0.001: (a) Fourier method; (b) Gaussian kernel method. 



reconstruction and the parameters using the mouse and the keyboard and to 
access to all output information of the method. 

The GUI has been realized in MATLAB (version 7.7.0 R2008b). Figure 
6.24 shows the main window of the GUI when it is started. Referring always 
to Figure 6.24 we can see it presents two windows where the original phantom 
and the reconstructed image of the phantom will be displayed. On the bot- 
tom of the window there is a panel that, thanks to pop- up menus, allows to 
choose the phantom and the reconstruction algorithm. There are three avail- 
able phantoms: the crescent-shaped phantom (introduced in section 2.7), the 
bull's eye phantom (Figure 6.25) and the Shepp-Logan phantom (Figure 2.5 
in section 2.5.1). Clicking on the options button a second window will be 
opened (hgura 6.26). Thanks to this options window it is possible to modify 
the number of sampled angle N, the number of samples on the t axis (2M+1) 
and the dimension of the output image (K 2 ). Moreover, depending on the 
selected algorithm, it is possible to change the predefined parameters used 
in the methods. For example, when using the back projection formula, one 
can choose both interpolation technique (nearest neighbor, linear, cubic) and 
the low pass filter (Ram-Lak, Shepp-Logan, cosine filter). Finally one can 
add a certain amount of noise to the Radon data to test the robustness of a 
method under the action of noise. Possible choices of noise are the Gaussian 
noise (with mean and variance decided by the user), Poisson noise and shot 
(or "salt and pepper") noise. Figure 6.27 shows the result of applying the 
back projection formula to the Shepp-Logan phantom adding Gaussian noise 
to data. 
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Figure 6.24: Graphical user interface 




Figure 6.25: Bull's eye phantom 
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Figure 6.26: Options window 
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Figure 6.27: Shepp-Logan phantom reconstruction with back-projection for- 
mula and Gaussian noise-data (mean fi — 0, variance a = 0.001). 
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It is also always possible to decide if to plot in a new figure the sinogram of 
the phantom, i.e. the sampled Radon transform used for the reconstruction. 
Figure 6.28 shows the case of Kaczmarz's method. In addition, depending 



I Figures -Figure 2 
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Figure 6.28: Sinogram and sparsisty of the matrix system of the Kaczmarz 
method. 

on the algorithm a further plot will be displayed: 

• Using the back projection formula, the interpolation of the convolution 
between the low pass filter and the Radon transform is shown; 

• With the Kaczmarz's method, the sparsity of the system matrix (figura 
6.28); 

• In the case of kernel methods, the image of the sysem matrix. 

The reconstruction begins pushing the start button. Figures 6.29(a) and 
6.29(b) show two different applications of the GUI. During the computation 
of the solution, in the MATLAB command window, it is displayed the status 
of the process, e.g. for the Kaczmarz's method 

Start reconstruction. . 
Radon transform computed. . 
computing ART system. . 
applying Kaczmarz's method.. 
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Figure 6.29: Example of reconstruction using the graphical user interface, (a) 
Multiquadric reconstruction of the crescent-shaped phantom; (b) Kaczmarz's 
reconstruction of the bull's eye phantom. 



residual : 
0.5895 



done . 

Where residual indicates the norm of the difference Axk — b of the solution 
at the final iteration k. 

Closing the GUI window or saving the workspace from the menu, one can 
find all information about the solution and the algorithm used in the output 
structure out. This structure contains the following fields: 

• radon: sampled Radon transform of the phantom (N x (2M + 1) ma- 
trix) ; 

• reconstruction: Reconstructed image of the phantom (K x K matrix 

); 

• phantom: name of the used phantom (string); 

• algorithm: name of the used algorithm (string); 

• options: options used in the reconstruction (structure depending on 
the algorithm). 

The options structure contains information about the sampling (N, M, K), 
on the noise (type of noise and mean and variance in the Gaussian case) and 
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a logical value that determines if the sinogram is plotted or not. Additional 
parameters depend on the algorithm: 

• Back projection formula: interpolation technique and the low pass fil- 
ter; 

• Kaczmarz's method: relaxed parameter A and a maximum number of 
iterations and a tolerance that decide when to stop the reconstruction 
process; 

• Kernel methods: the name of the kernel used (Gaussian, multiquadric, 
inverse multiquadric), the shape parameter of the kernel and which the 
window function has been used. 

The graphical user interface is designed in a way that is easy to extend 
its functionalities for example adding new phantoms or new reconstruction 
methods. Another improvement that can be added is the possibility to di- 
rectly compare the results of a reconstruction using two different methods or 
different parameters. 



Conclusions 



In this thesis we studied the problem of clinical image reconstruction. The 
problem was faced both from an analytical point of view, considering the 
mathematical aspect of the problem and the classical methods used to solve 
it, and also with a numerical approach, implementing new algorithms to solve 
it and comparing the behavior of the different methods. 

In the hrst part we focused on classical Fourier based methods. These 
methods are founded on the back projection formula and its discretization . 
We saw that in a discrete context it is possible to obtain only an approximated 
solution because of the presence of noise and the constrain to have only a 
finite amount of data. 

In the second part we introduced a different approach for solving the 
image reconstruction problem, called ART. With this kind of methods the 
solution is obtained solving a linear system. In particular positive definite 
kernel can be used to this aim, provided a regularization of the Radon trans- 
form functional. 

The regularization we used is to multiply the kernel function by a window 
function so that the Radon transform of the product function is finite. Then, 
we realized these algorithms using particular kernel and window functions 
and studied their behavior in function of shape parameters. 

In the last part of the thesis we compared kernel based with Fourier based 
methods. We saw that the quality of the reconstruction of the two methods 
is similar, also in the case of noisy data. The main limits of kernel based 
methods are the computational time, that grows exponentially with the size 
of the problem, and a bound for the number of data usable, indeed the linear 
system involved in the problem can become huge. 

Possible improvement and further works consist in using other kind of 
kernels and window functions. In this case the main difficulties can be finding 
the analytical expression of the Radon transform, or try other regularization 
techniques for the Radon transform integral. Implementing faster algorithms 
for solving the linear system, for example generating structured or sparse 
matrix, can be also another improvement. Moreover, an accurate study of the 
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approximation error can give useful informations, e.g. in the determination 
of optimal shape parameters. Finally, one can introduce a polynomial term 
in the expression of the solution and then use conditionally positive definite 
kernels. 

The big number of applications and the vastness of possibilities that can 
be followed in using the kernel based approach show why this research field 
has become so important in the last years. 



Appendix A: Inverse 
multiquadrics kernel matrix 



We compute the elements of the matrix A of the inverse multiquadrics re- 
construction problem (section 5.2). We recall that in that case 



with 

Ci = emax (-L, — \a\\/H 2 — r 2 + b) c 2 = emin (L, \a\ V H 2 — r 2 + b). 
Hence, all we have to do is the compute 



/ asinh (vT 



du, M > 0. 



+ u 2 

Using the logarithmic representation of asinh we can write 

/ = J log (VM 2 - u 2 + Vl + M 2 ) du-^J log (1 + u 2 ) du = 



— h — (^atanu + — log (1 + u ) — uj . 
Integrating Ii by parts 

log(VM 2 " — u 2 + VT+ M 2 ) + 

+ i (M 2 - u 2 ) - VM 2 - vt^/NPTl 
= ulog (VM 2 - u 2 + Vl + M 2 ) + I 2 . 
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Adding and subtracting M 2 in the numerator, we get 

M 2 -u 2 



I 



(M 2 - u 2 ) - VM 2 - u 2 VM 2 + 1 



du+ 



+ M 2 [ ; 1 du = 

J (M 2 - u 2 ) - VM 2 - u 2 VM 2 + 1 

= h + M 2 h. 

Setting a = acos-^ and c = \/l + m hi we have 

f a/M 2 -u 2 , , , f sin 2 a , 

h = ~ / / / du = Ml da = 

J VM 2 -u 2 + VM 2 + 1 J sina + c 

da 

+ c 



, , /• sin a - cr , , , 9 /• 1 

M / da + Mc 2 / 

j sm a + c 7 sin a 

2atan f , ^3 ) + a 

; \sina+(vc— l+vc+l) 2 / 



M(-cosa - ca) + Mc 



-u - (VM 2 + l) acos ^— J + 



2atan . = . + acos . 

VM 2 -u 2 + VM 2 + 1 + 1 J ^MJ 



u \ 



+ (M 2 + 1) 
Finally I 4 : 

I 4 = J — --<Jit - atan | si\:-r-:. I - ataiu/ 



1 , . / M 2 + l 

. , du = atan u\ — — 

(M 2 - u 2 ) - VM 2 - u 2 ^M 2 + 1 \ V M 2 - u 2 



Putting together the results, since I — Ii — I 3 — M 1" 4 , we obtain 



f , / /M 2 -u 2 \ 

i asmh ^VTT^J 



c?u = -asinh ( W — ) — (1 + M 2 )atanu+ 

2 \ V 1 + u 1 



+ JW+~1 (VJP + l - 1) acos (£) + M 2 atan ( u J ^£±1 ) + 



+ 2(M 2 + l)atan f . _ n ] . 

V 7 Vv / M 2 -n 2 + v / M 2 + l + l/ 

Where, of course, this formula is valid for \u\ < M. 



Appendix B: Compactly 
supported kernel matrix 



We compute the elements of the matrix A of the compactly supported func- 
tion reconstruction problem (section 5.4.1). We recall that in that case 



where 



a k ,j = / b(x s )(l - z/ 2 ||:r s ||) + (fs, 

JR 



b(x) 
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if t ^ 
if t = 0. 



Since ||x s || — r + s , one obtains 



a-ki = (1 — v 2 r 2 ) i b(x s )ds — u 2 / b(x s )s 2 ds 



(1 - v z r 2 )h - v 2 h if \r\ < - 

if Irl > - 

v 
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Setting t — x s ■ v = as + b and D s = [s : \s\ < \J ^ — r 2 , e\as + b\ < 1). 
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^- (2u 2 + l) 



(in /" C2 2 2 / 1 \ du 

/ -u acosh : — - — , 

ea Jet e \ ii / ea 
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where u = e(as + b) and 



Ci = max — 1 , be — e \ a 



r 2 , c 2 = min 1, be + e\a\ \ — — r 2 



thus 



6e 2 a 



e 2 a 



[3 arcsin u + u\/l — u 2 (2u 2 + 1)] c 2 + 



1 u 3 . ( 1 

- arcsin u H acosh - — . 

6 3 \\u\ 



6 



1 

3Ai 



Cl 

1^ 3 

- arcsin u + uVl — u 2 (u 2 + -) — 2n 3 acosh 



and the second integral becomes 



= / 

J c 

-I 



C2 2 

ci e 

C 2 2 



Vl — u 2 . 9 ./u — eb 

^— 2n 2 + l 

3 \ ea 



(in 



ea 



-+ 



2 



ti — e&\ 
ea / ea 



3e 4 a 3 I 12 

,3 



— ( b 2 e 2 + — I arcsin u + -be(l — u 2 ) 2 + 
.2 V 10/ 9 V 7 



^ / 1 \ 

(6u 2 - Ihbeu + 10fcV)acosh ( — + 

30 \M/ 

VI - « 2 /20 



60 



19, 



16fen 4 + (106 2 e 2 + — )u 3 + 



U , 2 /,r.2 2 In 28, 

— tett 2 + (156 2 e 2 - -)u - — 



C2 



Finally, if z/|r < 1, we have: 



3e 2 a 



- arcsin nil — v 2 r 2 



(•1) 



+ v / T^ 2 ^(n 2 + ^)( Z / 2 r 2 )- 



+u 3 acosh 



u 



y 2 qi{u) 

5e 2 a 2 



u 2 q 2 (u) 4u 2 be 
10e 2 a 2 ~ 3e 2 a 2 

\ -| c 2 

2„2< 



(l-^ 2 ))+ (-2) 



2(1 - ^V) 



(.3) 
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APPENDIX B: COMPACTLY SUPPORTED KERNEL MATRIX 



where 

q x (u) = Qu 2 - 15bsu + 10b 2 e 2 

q 2 (u) = ^n 5 - IQbeu 4 + (I0b 2 e 2 + ^)u 3 + 

_ 14 6ra 2 (156 2 £ 2 _ 1) _ 28 fe ^ 
3 v T 3 

We observe that if < 1, then Ci < c 2 always holds. 
While for v\r\ > 1, a^- = 0. 

At last we observe that because of the term acosh(|u| _1 ) in (.3), we have 
to consider apart the cases C\ = 0, c 2 = 0. In these cases, it is easy to see 
that, because of continuity, it is sufficient to consider the limit of (.3) for 
u — > 0, so that u 3 acosh(|u| _1 ) — > 0. 
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